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Abstract 

An  approximately  globally  convergent  numerical  method  for  a  ID  coefficient 
inverse  problem  for  a  hyperbolic  PDE  is  applied  to  image  dielectric  constants 
of  targets  from  blind  experimental  data.  The  data  were  collected  in  the  field  by 
the  Forward  Looking  Radar  of  the  US  Army  Research  Laboratory.  A  posteriori 
analysis  has  revealed  that  computed  and  tabulated  values  of  dielectric  constants 
are  in  good  agreement.  Convergence  analysis  is  presented. 

(Some  figures  may  appear  in  colour  only  in  the  online  journal) 


1.  Introduction 

In  this  paper,  we  test  the  ID  version  [30]  of  the  numerical  method  of  recent  publications 
[5-12,  26,  27,  31,  47]  for  the  case  when  the  time-resolved  backscattering  electric  signal  is 
measured  experimentally  in  the  field.  Measurements  were  performed  by  the  Forward  Looking 
Radar  built  in  the  US  Army  Research  Laboratory  (ARL).  All  kinds  of  clutters  were  present  at 

6  Formerly  with  University  of  North  Carolina,  Charlotte,  NC,  USA. 

7  Author  to  whom  any  correspondence  should  be  addressed. 
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the  site  of  data  collection.  The  data  are  severely  limited.  The  goal  of  this  radar  is  to  detect  and 
possibly  identify  shallow  explosive-like  targets.  Prior  to  this  effort,  the  focus  of  the  ARL  team 
was  on  the  image  processing  rather  than  on  the  target  detection  and  identification  [36].  The 
current  data  processing  procedure  of  ARL  delivers  only  the  energy  information.  The  algorithm 
of  this  paper  computes  values  of  dielectric  constants  of  targets  using  those  data.  These  values 
represent  a  new,  surprising  and  quite  useful  dimension  of  information  for  the  ARL  team.  A 
hope  is  that  these  values  might  be  helpful  in  the  target  detection  and  identification  process. 

The  UNCC/ChalmersGU  team  has  worked  only  with  the  most  challenging  case 
of  blind  experimental  data.  ‘Blind’  means  that  first  computations  were  made  by  the 
UNCC/ChalmersGU  team  without  any  knowledge  of  correct  answers.  Next,  computational 
results  were  sent  to  the  ARL  team.  The  ARL  team  has  compared  a  posteriori  those  results 
with  the  reality  and  then  revealed  correct  answers  to  the  UNCC/ChalmersGU  team. 

In  the  above  cited  works,  a  new  numerical  method  was  developed  for  some 
multidimensional  coefficient  inverse  problems  (MCIPs)  for  a  hyperbolic  PDE  with  single 
measurement  data.  ‘Single  measurement’  means  that  either  only  a  single  position  of  the  point 
source  or  only  a  single  direction  of  the  incident  plane  wave  is  considered.  Because  of  many 
dangers  on  the  battlefield,  the  single  measurement  arrangement  is  the  most  suitable  one  for 
military  applications.  There  were  two  goals  of  those  publications. 

Goal  1.  To  develop  such  a  numerical  method,  which  would  have  a  rigorous  guarantee 
obtaining  a  good  approximation  for  the  exact  solution  of  a  coefficient  inverse  problem  (CIP) 
without  using  advanced  knowledge  of  either  a  small  neighborhood  of  that  solution  or  of  the 
background  medium  in  the  domain  of  interest. 

Goal  2.  This  method  should  demonstrate  a  good  performance  on  both  computationally 
simulated  and  experimental  data. 

It  is  well  known  that  it  is  enormously  challenging  to  achieve  both  goals  1  and  2 
simultaneously.  Three  substantial  obstacles  are  combined  here:  the  minimal  information 
content  due  to  the  single  source  only,  nonlinearity  and  ill-posedness.  Therefore,  it  was 
inevitable  in  the  above-cited  publications  to  introduce  some  natural  approximations.  Although 
those  approximations  cannot  be  rigorously  justified  sometimes,  still  they  have  resulted  in  the 
simultaneous  achievement  of  both  goals;  see  more  details  in  sections  4  and  6.  The  numerical 
methods  of  [5-12,  26,  27,  30,  31,  47]  use  the  structure  of  the  underlying  PDE  operator  rather 
than  a  least-squares  functional.  The  above  thoughts  are  reflected  in  the  following  statement 
of  the  review  paper  [12]  ‘The  key  philosophical  focus  of  our  review  is  the  above  point  about 
natural  assumptions/approximations  which  make  the  technique  numerically  efficient.’ 

Because  of  the  measurement  scheme  of  the  Forward  Looking  Radar,  the 
UNCC/ChalmersGU  team  had  only  one  time-resolved  curve  for  each  target.  This  curve  was 
the  voltage  averaged  over  some  readings  (section  7).  Only  one  component  of  the  electric  field 
was  generated  and  measured  in  the  backscattering  regime.  The  reality  is  3D,  and  the  electric 
field  propagation  is  governed  by  the  full  Maxwell  system.  However,  the  above  data  structure 
has  left  us  with  no  choice  but  to  model  the  process  by  a  ID  CIP  using  only  one  hyperbolic 
PDE. 

The  main  challenge  of  working  with  these  data  was  a  huge  misfit  between  them  and 
computationally  simulated  data,  e.g.  compare  figure  3(b)  with  figures  4(b),  (d)  and  (f). 
Namely,  the  experimentally  measured  curves  are  highly  oscillatory,  unlike  the  computed 
ones.  Therefore,  to  make  the  experimental  data  suitable  for  the  inversion,  we  apply  a  new 
data  pre-processing  procedure,  which  is  a  crucial  step.  This  procedure  was  unbiased,  since 
we  have  worked  with  the  blind  data  only.  See  [10,  26]  and  chapter  5  of  [6]  for  a  similar  data 
pre-processing  procedure  for  the  case  of  transmitted  experimental  data.  Both  these  procedures 
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are  based  on  the  intuition  only.  The  only  justification  for  both  is  the  accuracy  of  obtained 
results;  see  also  comments  in  section  9  below. 

There  are  some  classical  numerical  methods  of  solving  ID  CIPs  for  hyperbolic  equations; 
see,  e.g.  [19,  25,  35]  and  references  therein.  However,  because  of  many  uncertainties  in  our 
experimental  data,  including  the  above  misfit  and  the  3D  reality  versus  the  ID  mathematical 
model  (also,  see  subsection  7.2),  it  is  yet  unclear  how  those  methods  would  perform  for  our 
data.  This  question  deserves  a  separate  research  effort. 

There  are  also  MClPs  with  multiple  measurement  data.  ‘Multiple  measurements’  means 
that  the  data  are  generated  either  by  the  point  source  running  along  a  manifold,  or  by  the  incident 
plane  wave,  whose  direction  varies  within  a  certain  cone.  These  MClPs  have  applications  in 
medical  imaging.  We  refer  to  [1,  13,  16,  17,  21,  24,  37-40]  for  some  non-local  reconstruction 
techniques  for  these  MClPs.  An  analogue  of  the  technique  of  [5]  for  the  case  of  a  2D  CIP 
for  an  elliptic  PDE  with  the  source  running  along  a  straight  line  was  developed  [28,  41, 
42].  Validation  of  the  latter  for  the  case  of  experimental  data  with  the  application  to  stroke 
detection  in  brains  of  small  animals  was  done  in  [43].  We  point  out  that  one  of  the  main  keys 
to  the  success  of  numerical  results  of  [1,  16,  17]  for  the  reconstruction  algorithms  of  Novikov 
[37-40]  is  the  use  of  approximate  mathematical  models.  The  same  can  be  said  about  [28, 
41-43].  This  concept  is  similar  to  the  one  above. 

In  section  2,  we  pose  a  CIP  for  a  ID  wave-like  PDE.  In  section  3,  we  study  some 
properties  of  the  Laplace  transform  of  the  solution  of  the  forward  problem.  In  section  4, 
we  discuss  the  concept  of  the  approximate  global  convergence  property.  In  section  5,  we 
present  our  numerical  method.  Convergence  analysis  can  be  found  in  section  6.  In  section  7, 
we  describe  the  experimental  setup,  main  uncertainties  in  the  experimental  data  and  the  data 
pre-processing  procedure.  Imaging  results  are  presented  in  section  8.  Section  9  is  devoted  to 
discussion. 

2.  Statements  of  forward  and  inverse  problems 

Let  the  function  er(x),  x  e  R  be  the  dielectric  constant,  spatially  distributed  on  the  real  line. 


Let  the  number  d  >  1 .  We  assume  below  that 

er{x)  e  [1,  d),  x  e  [0,  1],  t,  e  C1  (M),  (2.1) 

e,.(v)  =1,  *  i  (0,  1).  (2.2) 

Thus,  the  interval  (0,  1)  is  our  domain  of  interest  in  the  inverse  problem.  The  forward 
problem  is 

er(x)utl  —  uxx,  rel,  t  e  ( 0,  oo),  (2.3) 

u  (.v,  0)  =  0,  U[  ( x ,  0 )  —  8  (x  —  a'o ) ,  (2.4) 

xq  —  const.  <  0.  (2.5) 

CIP.  Determine  the  coefficients  r{x ),  assuming  that  the  following  function  g(t)  is  known 
u(0,t)  =  g(t),  te(  0,  oo).  (2.6) 


The  function  g(t)  models  the  backscattering  data  measured  by  the  Forward  Looking 
Radar.  The  condition  re  (0,  oo)  in  (2.6)  is  not  a  restrictive  one  from  the  numerical  standpoint. 
Indeed,  we  use  the  Laplace  transform  (2.7)  to  solve  this  CIP.  Since  the  kernel  e-st  of  this 
transform  decays  rapidly  with  respect  to  r,  then  in  actuality  a  finite  time  interval  is  sufficient. 
In  addition,  the  data  resulting  from  the  data  pre-processing  procedure  have  a  finite  support  in 
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(0,  oo) ,  see  figures  5(a)  and  (b).  Since  this  is  the  ID  case,  uniqueness  theorems  for  this  CIP 
are  well  known;  see  e.g.  [35]  and  references  therein. 

Consider  the  Laplace  transform 

poo 

w(x ,  s)  —  u  (x,  t)  e~stdt  :=  Cu,  s  ^  s  =  const.  >  0,  (2.7) 

Jo 

where  the  number  s  —  s  (er)  is  sufficiently  large.  We  call  the  parameter  s,  pseudo  frequency. 
Hence, 

u>xx  —  s2er(x)w  =  —  S  (x  —  xo),  x  e  R,  (2.8) 

lim  w(x,  s)  =  0;  (2.9) 

|a:| — >oo 

see  theorem  3.1  for  the  proof  of  (2.9).  By  (2.6) 

w  (0,  s)  =  (pis )  :=  Cg,  s  ft  s.  (2.10) 

We  also  need  to  know  the  derivative  wx  (0,  ,v), 
wx  (0,  s)  =  p(s),  s  ^  s. 

Consider  the  fundamental  solution  w o(x,  ,v  )  of  problem  (2.8)  and  (2.9)  for  e,-(x)  =  1.  Then 


u’o(x,s)=  (2j)_1  exp  (— s  \x  —  xo|).  (2.11) 

Denote  w(x,  i)  =  w(x,  s )  —  wq(x ,  s).  Using  (2. 8)— (2. 10),  we  obtain 

wxx  —  s2e,fx)w  =  s2  (e,  (x)  —  1)  wq,  x  e  K,  (2.12) 

lim  ui(x,  ^)  =  0,  (2.13) 

|x| — >CxD 

w  (0,  s )  =  <p(s)  :=  (p(s)  —  (2i)_1  exp  (— s  |xo|) .  (2.14) 

Also,  (2.2)  and  (2. 12)-(2. 14)  imply  for  x  <  0, 

Wxx  —  s2w  =  0,  x  <  0,  (2. 15) 

lim  w(x,s)  —  0,  (2.16) 

x—> — OO 

w  (0,  .S’)  =  ip(s).  (2.17) 

It  follows  from  (2. 1 5)— (2. 1 7)  that  u>(x,s )  =  <p(s)esx,  x  <  0.  Hence,  wx  (0,  .s)  =  sip{s )  = 
s((p(s)  —  (2s)_1  exp  (—s  |xo|))  and  thus  since  xo  <  0,  we  have 

wx  (0,  s)  :=  p{s)  =  scpis)  —  exp(.?xo).  (2.18) 


3.  Some  properties  of  the  function  wix,  s) 

Below  C“,  a  e  (0,  1)  are  Holder  spaces.  In  this  section,  we  establish  some  properties  of  the 
function  w  =  Cu,  which  we  need  for  the  convergence  analysis.  Let  f  it),  t  >  0  be  a  piecewise 
continuous  function  such  that  |  /  (r)  <  Cza',t>  1,  where  C  =  C  if),  a  —  a  if)  =  const.  > 
0.  Consider  two  types  of  Laplace  transforms, 

£[(/)w=v^i 

poo 

Cl  if)  is)  =  fit)  e~A  dr,  s>  y/aif). 

Jo 
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Lemma  3.1.  ZL  [£1  (/)]  (j)  =  C  (/)  ( s )  for  s  >  a  (/),  where  the  operator  C  was  defined 
in  (2.7). 


Proof.  Using  formula  (28)  of  section  4.5  of  the  book  [2],  we  obtain 
1 


n 


VJFf3/2 


r  exp 


dr  =  e~ST ,  'is  >  0. 


□ 


Theorem  3.1.  Consider  the  function  er(x )  e  C“  (R)  satisfying  the  rest  of  conditions  (2.1) 
and  (2.2).  Let  (2.5)  hold.  Then  for  any  number  s  >  0  there  exists  a  unique  solution 
p(x,  s)  e  C2+a  (|x  —  xq|  ^  f)  fl  C  (R) ,  V/3  >  0  of  the  following  problem: 


Pxx  —  S2Sr(x)p  =  —8  (x  —  X'o),  x  e  R, 

(3.1) 

lim  p(x,  i)  =  0. 

(3.2) 

|jc| — >-oo 

0  <  p(x,  i)  ^  wq(x ,  i),  Vx  e  R. 

(3.3) 

In  addition,  let\x  —  xol  7s  f  —  const.  >  0.  Then  there  exists  a  sufficiently  large  number 
T=~s~(d,  ft)  such  that 

p(x,  s)  >  Wd(x,  s),  Vx  e  [xo  +  f,  00),  Vi  ^  s',  (3.4) 

where  Wd(x ,  s)  is  the  fundamental  solution  of  equation  (2.8)  for  the  case  er(x )  =  d , 

Wd(x,s)  =  (2s)_1  exp(— sVd\x  —  xol)-  (3.5) 

Finally,  there  exists  a  sufficiently  large  number  s  —  s  (er)  >  0  such  that 

w(x,  s)  C  (u)  (x,  i)  =  p(x,  s),  Vi  ^  s  (er) ,  Vx  e  M.  (3.6) 

Thus,  (3.2)  implies  that  (2.9)  holds  for  s  7?  s  =  s  (er). 


Proof.  Let  the  function  v  (x,  t)  be  the  solution  of  the  following  Cauchy  problem: 


er(x)v,  =  vxx, 

(x,  r)  e  R  x  (0,  oo), 

(3.7) 

v  (x,  0)  =  8  (x  - 

-x0). 

(3.8) 

Let  the  function  Uq  (x,  t)  be  the  solution  of  problem  (3.7),  (3.8)  for  the  case  er(x)  =  1, 


v0  (x,  r) 


1 


:  exp 


(x  -  xo£ 

4r 


2Vrrr 

Let 7)  (x,  t)  —  v  (x,  t)  —  vq  (x,  r).  Then 

vxx  ~  e,  (x)Vt  =  (£r(x)  -  1)  v0t,  v  (x,  0)  =  0. 


(3.9) 


Detailed  estimates  of  the  fundamental  solution  of  a  general  parabolic  equation  in  [32,  chapter 
4]  imply  thatu  e  c2+“’1+“/2  (M  x  [0,  T]),  WT  >  0.  Denote 


v  (x,  t)  =  /  fi  (x,  r)  dr  =  / 

Jo  Jo 


(v  —  Vq)  (x,  r)  dr. 


By  (3.9) 


Vxx  —  Sr(x)v,  =  (£r(x)  —  1)  Do,  v  (x,  0)  =  0. 


(3.10) 


(3.11) 
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By  (2.2)  and  (2.5),  er(x)  —  1  =  0  in  a  neighborhood  of  the  source  point  { a'q } .  Hence, 
applying  (2.2),  we  obtain  (er(x)  —  1)  i>o  (x,  t )  ^  0  in  K  x  [0,  oo) .  Hence,  applying  the 
maximum  principle  of  [23,  theorem  1  of  chapter  2]  to  (3. 1 1)  we  obtain 


v  (x,  t)  ^  0  in  lx  [0,  oo). 


(3.12) 


On  the  other  hand,  theorem  1 1  of  chapter  2  of  [23]  ensures  that  the  fundamental  solution  of 
the  parabolic  equation  is  positive  for  t  >  0.  Hence, 


v  (x,  t)  >  0,  t  >  0. 

Hence,  (3.10),  (3.12)  and  (3.13)  imply  that 


0  < 


I  v  (x,  r)  dr  ^  1  Vo  (x,  r)  dr. 

Jo  Jo 


(3.13) 


(3.14) 


Using  the  Fubini  theorem,  formula  (27)  of  section  4.5  of  the  book  [2]  and  (2. 1 1),  we  obtain 


r  (  ('  ,  ^  a\  1  ^  ,  ,  m(x,s) 

£2  I  J  V0  (X,  r)  dr  1  =  -j£2  (v0)  =  — ^ — • 

Hence,  using  (3. 13)  and  (3. 14)  and  Fubini  theorem,  we  obtain 


(3.15) 


£2 


[‘  ,  ,  \  1  .  ,  ,  ,  ,  y(x,s)  w0{x,  s) 

/  v  (x,  r)  dr  1  =  —  £2  (v)  (x,  s)  :=  — <  - y — , 

Jo  J  s  s  s 


(3.16) 


By  (2.11),  (3. 13)  and  (3.16) 

lim  y(x,  s )  =  0. 

|jc| — >00 

Next,  by  (3.10)  and  (3.11)  vxx  =  sr(x)v  —  Vo ■  Hence,  | T7^|  <  dv  +  vq.  Therefore, 

nOO 

/  \Vxx  (x,  f)|  e-'2'  d t  <  00. 

Jo 

Hence, 


(3.17) 


3; 


_ 2  /'00~ 

/  v  (x,  t)  e  st  d t  =  I  vxx  (x,  t)  e  s~' d t  =  sr(x)y  —  wq. 
Jo  Jo 


On  the  other  hand, 

vxx  O,  t)  e 


L 


vxx  (x,  r  )  dr 


dr 


fd 


Uoxc  (x,  r)  dr 


(3.18) 


dr 


=  s  2  (yxx-  Woxx)  ■ 

Comparing  this  with  (3. 18)  and  keeping  in  mind  that  Wqxx  —  s2wq  =  —  8  (x  —  xo) ,  we  obtain 

yxx  ~  s2sr(x)y  =  -S  (x  -  x0).  (3. 19) 

By  (3. 17)  and  (3. 19),  the  function  y(x,  s )  =  £2  (v)  (x,  s)  satisfies  conditions  (3.1)  and  (3.2).  It 
follows  from  above  that  y(x,  s)  e  C2+a  (|x  —  xo|  ^  /3)  DC  (K),  Vj8,  s  >  0.  Uniqueness  of  the 
solution  of  problem  (3.1)  and  (3.2)  for  this  class  of  functions  can  be  easily  proven  in  the  standard 
way  using  the  maximum  principle  for  elliptic  equations  and  (3.2).  Hence,  problem  (3.1) 
and  (3.2)  have  the  unique  solution  p(x,  s )  :=  y(x,  s )  e  C2+a  (|x  —  xo|  ^  jS)DC  (M),  Vj3,  s  >  0. 
The  right  estimate  (3.3)  follows  from  (3.16),  and  the  left  estimate  follows  from  (3.16). 
Lemma  3.1  implies  (3.6). 

We  now  prove  (3.4).  We  obtain  from  (2.2),  (3.1)  and  (3.2): 
exp  (— s  lx  —  Xnl)  s 

p(x,  s )  =  ' - —  -  -  /  exp  (-S  \x  -  §1)  (er  (?)  -  1  )p  (?,  s)  d?.  (3.20) 

2s  2  Jo 
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Hence,  by  (2.1)  and  the  right  inequality  (3.3) 
exp  (— s  |x  —  xol)  (d  —  1) 


P(X-S)>  2,  4 

Assume  now  that  x  e  (xq,  0).  Then, 


f 


exp  (—s  |x  —  f  |)  exp  (—s  |£  -x0|)  d£. 


p(x,  s)  ft 


exp[— s(x  —  x0 ) ]  ( d  —  1)  exp  [j  (x  +  x0)] 

2s  4 

exp  [— s(x  —  xq)] 


/' 


,-24 


2s 


\  _  (d  —  1 )  ( 1  —  e~2;) 
4 


Jo 

a  2  SX 


Therefore, 


P("]'  S\  >  exp[i(v7  -  l)(x  -  x0)] 
wd(x,  s) 


1  - 


(d-  1)(1  -e~2s)  e2“ 


Vx  e  (xq,  0). 


Hence,  if  ft  =  const,  e  (0,  |xol),  then  there  exists  a  number  T=  sill,  ft)  such  that 
p(x,  s) 


>1,  Vi  ^  s  =  s  (d,  ft) ,  Vx  e 


xo  +  -.-*o  +  , 


wd(x,  s) 

Consider  now  the  difference  p(x,  s)  =  p(x,  s)  —  wd(x,  s).  Then, 

Pxx  —  S2sr(x)p  =  s2  (s,  (x)  —  d)\Vd  <  0,  x  e  [xq  +  ft,  oo), 


lim  p(x,  s)  =  0. 

x—>oo 


(3.21) 

(3.22) 

(3.23) 


Let  a  >  xo  +  ft  be  an  arbitrary  number.  Consider  the  function  p(x,  s)  on  the  interval  [xo  +  ft,  a] 
fori  Js  's  (d,  ft).  Then,  the  maximum  principle  for  elliptic  equations  (3.21)  and  (3.22)  implies 
that  the  negative  minimum  of  p(x,  s)  on  this  interval  can  be  achieved  only  at  x  —  a  Setting 
a  ->  oo  and  using  (3.23),  we  obtain  that  p(x,  s)  ft  0  for  x  e  [xo  +  ft,  oo),  s  ft's  (d,  ft). 

Assume  now  that  there  exist  a  point  x  e  [xo  +  ft,  oo)  and  a  number  i  ^  T(d,  ft)  such  that 
p  (x,  i)  =  0.  Since  p(x,  s)  ft  0  for  x  e  [xo  +  ft,  oo),  therefore, 

p  (x,  j)  =  0  =  min  p(x,  i).  (3.24) 

[.*o+£°°) 

Hence,  p^  (x,  s)  ft  0.  However,  this  and  (3.24)  contradict  to  the  inequality  in  (3.22).  □ 


Corollary  3.1.  Let  the  function  sr  (x)  e  Ca  (R)  satisfy  the  rest  of  conditions  (2.1)  and  (2.2).  For 
each  s  >  0  and  for  each  xel,  the  integral  of  the  Laplace  transform  £2  (v  (x,  t))  converges 
absolutely. 


Proof.  See  (3. 14),  (3. 15)  and  (3. 16). 


□ 


Remark  3.1.  By  (3.6)  the  function  ip(s)  in  (2.10)  is  defined  only  for  5  ft  s(er).  However, 
using  lemma  3.1  and  corollary  3.1,  we  set  below  ip(s)  :=  £2  [£1  (g)]  (s),  Vs  >  0,  where  the 
function  g  ( t )  is  defined  in  (2.6). 

Theorem  3.2.  Let  the  function  er(x)  e  Ca  (R)  satisfy  the  rest  of  conditions  (2.1),  (2.2) 
and  (2.5).  Let  the  number  s  ft  s\|xo|,  d),  where  the  number's  was  defined  in  theorem  3.1.  Let 
the  function  p  (x,  s;  er)  p  (x,  s)  e  C2+a  (|x  —  xol  ^  ft)  fl  C  (R),  V/3  >  0  be  the  solution  of 
the  problem  (3.1)  and  (3.2).  Denote 

fk  (x,  s\  er)  —  s~2 3*  [In  p  (x,  s)] ,  k  =  1,  2. 
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Then,  there  exists  a  constant  B  —  B  (.To,  d,  s)  >  1  such  that  for  all  such  functions  er 


II  fk  (,X,  S',  £r)llc[0,l] 


^ B ,  £=1,2. 


Also,  for  any  two  functions  and  ef2>  satisfying  conditions  of  this  theorem 

II/*  ( * -  £rV)  ~  fk  iX •  *  £r2’)  L(0.n  <B\\e> 


(1) 


,(2)| 


II  £2(0,1)’ 


£=1,2. 


(3.25) 


(3.26) 


Proof.  In  this  proof  B  =  B  (xq,  d,  s)  >  1  denotes  different  constants  depending  on  the  listed 
parameters.  We  have 


,  _  Px 

J 1  -2  ’ 

£  P 


f  Pxx  Px 
J2  _2  — 2  o 

£  P 


s2p2 


Hence,  by  (3.4) 


(3.27) 


l/l  (x,  s;  Er)  I  <  B  \px  (t,  5;  £r)|  , 


1/2  (x,  .S';  e,.)|  <  B[|p.«|  +  \px\2]  (x,  s\  er) ,  x  e  [0,  1]. 


(3.28) 


To  estimate  the  function  \px\,  we  use  the  integral  equation  (3.20).  We  have  for  t  e  (0,  1), 

Px(x,s)  =  _exP[  s(x — ^11  +  S—  /  sgn  (x  —  £)  exp  (—s  |t  —  §|)  (sr($)  -  l)p($,s)  d$. 

Hence,  (2.1)  and  (3.3)  imply  that  \px\  ^  B  for  x  e  [0,  1].  Since  S(x  —  To)  =  0  for 
t  e  [0,  1],  (3.3)  implies  that  pxx  =  s2sr(x)p  for  x  e  [0,  1].  Hence,  by  (2.1)  and  (3.3) 
|  Pxx  I  <fi,T£  [0,  1],  Thus,  (3.28)  implies  (3.25). 

We  now  prove  (3.26).  Let  p(x,  s)  =  p  (t,  s;  £*n)  —  p  (t,  5;  £^2))-  Then  by  (3.1) 

Pxx  —  s2ell)(x)p  =  s2  (e^  —  e®)  (x)p  (x,  5;  sj.2>) ,  ret.  (3.29) 

In  addition,  it  follows  from  (3.3)  and  (3.20)  that  functions  d]xp( x,  s) ,  j  =  0,  1  and  2  decay 
exponentially  with  |t|  —*■  oo.  Hence,  multiplying  (3. 29)  by  p,  integrating  over  R  and  using  (2. 1) 
and  (3.3),  we  obtain  ||p||Hi(R)  <  B  ||e'n  —  s^2>  ||  Q  u-  Hence,  (3.25)  and  (3.29)  lead  to 

IIpIWd^^I^-^Ud- 

Thus,  (3.4)  and  (3.27)  imply  (3.26).  □ 


4.  Approximate  global  convergence 

4.1.  The  concept 

As  to  the  ID  version  of  the  technique,  which  was  published  in  [30]  and  which  is  used  here, 
originally  the  work  of  [30]  was  considered  ‘only  as  a  preliminary  step  before  applying  similar 
ideas  to  2D  and  3D  cases’  [30,  p  125].  The  authors  of  [30]  meant  the  incorporation  of  the 
quasi-reversibility  method  (QRM)  in  the  techniques  of  [5-12,  26,  27,  31,  47];  see  section  5. 

Least-squares  functionals  for  MCIPs  with  single  measurement  data  suffer  from  multiple 
local  minima  and  ravines.  This  implies  the  local  convergence  of  conventional  numerical 
methods  for  those  MCIPs.  The  latter  mean  that  iterations  should  start  from  a  point  located  in 
a  sufficiently  small  neighborhood  of  the  exact  solution.  The  central  question  addressed  in  the 
publications  above  was:  How  does  one  construct  an  effective  numerical  method,  which  would 
lead  to  a  point  in  a  sufficiently  small  neighborhood  of  the  exact  solution  of  an  MCIP  without 
any  a  priori  information  about  that  neighborhood  with  a  rigorous  guarantee  of  reaching  that 
neighborhood ? 

Because  of  the  enormous  challenge  of  addressing  this  central  question  (section  1),  the 
approach  of  [5-12,  26,  27,  30,  31,  47],  as  well  as  that  of  the  current  paper,  consist  of  the 
following  six  steps. 
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Step  1.  A  reasonable  approximate  mathematical  model  is  proposed.  The  accuracy  of  this 
model  cannot  be  rigorously  estimated. 

Step  2.  A  numerical  method  is  developed,  which  works  within  the  framework  of  this  model. 
Step  3.  A  theorem  is  proven,  which  guarantees  that,  within  the  framework  of  this  model, 
the  numerical  method  of  step  2  indeed  reaches  a  sufficiently  small  neighborhood 
of  the  exact  solution,  as  long  as  the  errors  both  in  the  data  and  in  some  additional 
approximations  are  sufficiently  small.  It  is  central  to  our  approach  that  this  theorem 
should  not  rely  either  on  the  assumption  about  a  knowledge  of  any  point  in  a 
small  neighborhood  of  the  exact  solution  or  on  the  assumption  of  knowledge  of 
the  background  medium  inside  of  the  domain  of  interest. 

Step  4.  That  numerical  method  is  tested  on  computationally  simulated  data. 

Step  5.  (Optional).  The  numerical  method  of  step  2  is  tested  on  experimental  data.  To  have  a 
truly  unbiased  case,  blind  data  are  preferable.  This  step  is  optional  because  it  is  usually 
not  easy  to  actually  obtain  experimental  data. 

Step  6.  Finally,  if  the  results  of  step  4  and  (optionally)  step  5  are  good  ones,  then  that 
approximate  mathematical  model  is  proclaimed  as  a  valid  one. 

Step  6  is  logical  because  its  condition  is  that  the  resulting  numerical  method  is  proved  to  be 
effective,  i.e.  goal  2  is  successfully  met;  see  also  remark  4.4.  After  a  finite  (rather  than  infinite) 
number  of  iterations  it  is  sufficient  to  achieve  the  small  neighborhood  of  the  exact  solution. 
Next,  because  of  approximations  in  the  mathematical  model,  the  resulting  solution  can  be 
refined  via  a  locally  convergent  numerical  method.  The  latter  led  to  a  two-stage  numerical 
procedure  in  [6-12,  31,  47].  In  the  first  stage,  the  numerical  method  of  step  2  is  applied.  In  the 
second  stage,  the  adaptive  finite  element  method  (adaptivity)  for  MCIPs  takes  the  solution  of 
the  first  stage  as  the  starting  point  for  iterations  and  refines  it. 

4.2.  Definition 

Considerations  of  subsection  4. 1  led  to  the  introduction  of  the  approximate  global  convergence 
property  in  [6,  12].  Since  this  definition  is  new,  it  is  worth  repeating  it  here.  Recall  that  one  of 
the  backbone  principles  of  the  theory  of  ill-posed  problems  is  that,  given  an  ill-posed  problem, 
it  is  reasonable  to  assume  uniqueness  and  existence  of  the  exact  solution  of  this  problem  for 
the  ‘ideal’  noiseless  data  [6,  22,  46].  It  is  not  necessary  for  an  algorithm  addressing  the  above 
central  question  to  start  from  any  such  point.  In  fact,  it  is  sufficient  for  it  to  start  from  such 
a  reasonable  point,  which  would  not  contain  information  about  a  small  neighborhood  of  the 
exact  solution.  The  question  of  the  validity  of  the  approximate  mathematical  model  M  of 
definition  4. 1  should  be  addressed  in  steps  4  and  5  of  subsection  4. 1 . 

Definition  4.1.  (Approximate  global  convergence)  [6,  12],  Consider  a  nonlinear  ill- 
posed  problem  P.  Suppose  that  this  problem  has  a  unique  solution  x*  e  B  for  the 
noiseless  data  y*,  where  B  is  a  Banach  space  with  the  norm  ||-||B.  We  call  x*  the 
‘exact  solution  ’  or  the  ‘correct  solution  ’.  Suppose  that  a  certain  approximate  mathematical 
model  M  is  proposed  to  solve  the  problem  P  numerically.  Assume  that ,  within  the 
framework  of  the  model  M,  this  problem  has  a  unique  exact  solution  x*M  and  let 
x*M  =  x*.  Consider  an  iterative  numerical  method  for  solving  the  problem  P.  Suppose 
that  this  method  produces  a  sequence  of  points  [x„}^=1  C  B,  where  N  e  [l,oo).  Let 
a  sufficiently  small  number  9  e  (0,  1).  We  call  this  numerical  method  approximately 
globally  convergent  of  the  level  9,  or  shortly  globally  convergent  if  within  the  framework 
of  the  approximate  model  M ,  a  theorem  is  proven,  which  guarantees  that,  without  any 


9 


Inverse  Problems  28  (2012)  095007 


A  V  Kuzhuget  (UNCC/ChalmersGU  team)  et  at 


a  priori  knowledge  of  a  sufficiently  small  neighborhood  ofx* ,  there  exists  a  number  N  e  [  1 ,  A^) 
such  that 


\xn  —  x*\\B  <  0 ,  V«  >  TV. 


(4.1) 


Suppose  that  iterations  are  stopped  at  a  certain  number  k  ^  TV.  Then,  the  point  Xk  is  denoted 
Xk  :=  Xgiob  and  is  called  ‘the  approximate  solution  resulting  from  this  method’. 

Remark  4.1.  We  repeat  that  we  have  introduced  this  definition  because,  in  simple  terms, 
nothing  else  works  for  MCIPs  with  single  measurement  data;  see  comments  about  goals  1  and 
2  in  section  1 . 

Remark  4.2.  The  most  important  requirement  of  definition  4. 1  is  that  this  numerical  method 
should  provide  a  sufficiently  good  approximation  for  the  exact  solution  x*  without  any 
a  priori  knowledge  of  a  sufficiently  small  neighborhood  of  x*.  Furthermore,  one  should 
have  a  rigorous  guarantee  of  the  latter  within  the  framework  of  the  model  M.  In  other  words, 
step  1  of  subsection  4.1  should  be  addressed. 

Remark  4.3.  Unlike  the  classical  convergence,  this  definition  does  not  require  x„  =  x*. 

Furthermore,  the  total  number  TV  of  iterations  can  be  finite. 

Remark  4.4.  As  to  the  use  of  the  approximate  mathematical  model  M,  all  equations  of 
mathematical  physics  are  approximate  ones.  The  main  criterion  of  their  validity  is  the  accuracy 
of  descriptions  of  experimental  data,  which  is  exactly  what  we  do.  Also,  it  is  well  known  that 
the  Huygens-Fresnel  optics  is  not  yet  rigorously  derived  from  the  Maxwell  equations;  see 
section  8.1  of  the  classical  book  of  Born  and  Wolf  [14].  Nevertheless,  it  is  the  Huygens- 
Fresnel  theory  which  describes  the  experimental  data  of  the  diffraction  optics  very  well. 
Furthermore,  the  entire  optical  industry  nowadays  is  based  on  the  Huygens-Fresnel  theory. 
On  the  other  hand,  the  derivation  of  this  theory  from  the  Maxwell  equations  is  based  on  some 
non-rigorous  approximations.  Analogously,  although  the  numerical  method  of  [5-12,  26,  27, 
31,  47]  is  based  on  an  approximate  model,  goal  2  has  been  consistently  achieved. 

5.  Numerical  method 

Although  a  detailed  description  of  the  numerical  method  of  this  section  can  be  found  in  [30], 
we  still  present  it  rather  briefly  here  in  order  to  refer  to  some  formulas  in  the  convergence 
analysis  in  section  6.  Also,  because  of  (2.2),  we  consider  here  the  function  In  ( w/wq )  instead 
of  In  w(x,  s)  in  [30].  Let  w(x,  5)  e  C2+a  (|x  —  xq\  ^  f)C\  C  (R),  V/J,  s  >  0  be  the  solution  of 
problem  (2.8)  and  (2.9);  see  theorem  3.1. 

5.7.  Integral  differential  equation 

Since  by  (3.3)  w(x,  s)  >  0,  we  can  consider  the  function  r(x,  s), 


r(x,s)  =  s  2  [In  w(x,  s)  —  In  wq(x,  s)]  =  5  1 1n  [(w/wo)  (x,  s)] . 


Then,  (2.8),  (2.10)  and  (2.18)  imply  that 
rxx  +  s2r^  —  2 srx  =  er(x)  —  1,  x  >  0, 

r  (0,  s)  =  (fio(s),  rx  (0.  s)  =  cpt  (s), 


(5.1) 


(5.2) 


<p0(s)  =  s  2  [ln</)(s)  —  ln(2s)]  +  xqs  \  (fi\(s)  =  2/s  —  expCsxoHs'VU))  ’.  (5.3) 
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By  (2.2),  (2.5),  (2.8)  and  (2.9),  w(x,  5)  =  C(s)e  “for* ^  1,  where  C(s)  is  a  certain  function 
s  dependent  on  s.  Hence, 

r,(l,  ■*)  =  <).  (5.4) 

Differentiate  equation  (5. 1)  with  respect  to  5.  Then, 

q(x,  s)  =  dsr(x,  s),  1 AoCs)  =  <p'0(s),  1//1  (,s)  =  cp[{s),  (5.5) 


r(x,  s)  =  —  J  q  (x,  r)  dr  +  V  ( x ,  s), 

V  (x,  s)  =  s~2  [In  w  (x,  s )  —  In  wq  (x,  5)]  =  r  (x,  s). 


(5.6) 

(5.7) 


Here,  s  >  0  is  a  sufficiently  large  number  which  is  chosen  in  numerical  experiments.  We 
call  V  (x,  s)  the  tail  function.  Actually,  s  is  the  regularization  parameter  of  this  method. 
Using  (5. 1)— (5.7),  we  obtain 

-,2 


l 


qxx-2s~qx  /  qx  (x,  r)  dr  +  2s  qx  (x,  r)  dr 


u; 

i: 


i 


-  2 sqx  +  2  qx  (x,  r  )  dr 


+  2 s2qxVx  -  4 sVx  I  qx  (x,  r )  dr  +  2s  (Vv)2  -  2VX  =  0,  s  e  [5,  5],  (5.8) 


q  (0,  s)  =  Vfo(>y)>  Qx  (0,  .S’)  =  f\  is),  qx  (1,  s)  =  0,  se[s,4 


(5.9) 


Lemma  2.1  of  [30]  implies  for  k  —  0,  1  and  2 

1 


Dkxr(x,  s)  =  Dkx 


P 

Jxo 


Sr  (?)  d§  -  (X  -  X0) 


o  - 


x  >  0,  s  —*■  00.  (5. 10) 


Hence,  the  function  V  (x,  s)  =  r  (x,  .v)  is  small  for  large  values  of  s. 


l|V(x,5)||<?B,,1]  =  0(1/5).  (5.11) 

The  main  difficulty  of  this  numerical  method  is  in  the  solution  of  problem  (5.8)  and  (5.9). 
Equation  (5.8)  has  two  unknown  functions  q  and  V.  To  approximate  both  of  them,  we  use 
a  predictor/corrector-like  scheme.  First,  given  an  approximation  for  V,  we  update  q  via 
equation  (5.8)  and  (5.9).  Next,  we  update  the  unknown  coefficient  sr(x)  and  solve  the  forward 
problem  (2.8)  and  (2.9)  for  s  :=s  with  this  updated  coefficient  sr(x).  This  is  our  predictor-like 
step.  On  the  corrector-like  step  we  update  the  tail  function  V  (x,  s)  via  (5.7). 

Consider  a  partition  of  the  interval  [s,  5]  into  N  small  subintervals  with  the  grid  step  size 
h  >  0  and  assume  that  the  function  q  (x,  s)  is  piecewise  constant  with  respect  to  s, 


s  =  Sn<  sjv-  1  <  •  •  •  <  so  =  s,  Si- 1  —  Si  —  h\  q(x,  s)  =  qn(x),  for  j  e  ( sn ,  s„_i].  (5. 12) 

Let  fi  1  be  a  large  parameter  which  should  be  chosen  in  numerical  experiments.  Multiply 
both  sides  of  equation  (5.8)  by  the  Carleman  weight  function  (CWF), 

Cn,ii(.s)  =  [x(sn—  1  s)]>  S  G  (xn,  Sn_l],  (5.13) 

and  integrate  with  respect  to  j  e  (sn ,  ,y„_i ).  The  CWF  is  introduced  in  order  to  mitigate  the 
influence  of  the  nonlinear  term  Bn  (/i.  h)  (q^f  in  the  resulting  equation.  If  h  is  fixed,  then 

B„{ix.  h)  =  D(/x_1),  h  ->  00.  (5.14) 


We  ignore  the  nonlinear  term,  since  we  observed  in  our  computations  that  it  provides  only 
an  insignificant  impact  in  results  for  ji  f  50.  For  each  n,  we  perform  inner  iterations  with 
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respect  to  the  tail  function.  This  way  we  obtain  functions  qnk  and  Vn^.  The  equation  for  the 
pair  (q„'k,  Vn,k)  is 


n—  1 


n—  1 


n- 1 


%.k 


Ai.nh  4  -  Ai.nK,k  -  2A2,«  q'n,k  =  -A2,„h2  +2 4j 


3= 0 


J= 0 


j=o 


n — 1 


+  2A2 t„V^k  j  h  Y,  q'A  —  A2,„  «,)  +  2A2i„Vn'jjt,  go  :=  0,  xe  (0,  1), 

(5.15) 


<3n,k  (0)  =  fo,n,  q'n.k  (°)  =  ,»>  4n,k  C1)  =  0,  (5.16) 

\  r-'  1 

•n  =  yiJ  fo(s)ds,  fi  ,n=jiJ  ifi(s)ds.  (5.17) 

Here,  A\Jt  and  A2  „  are  certain  numbers,  whose  exact  expression  is  given  in  [5,  6].  It  is  known 
that 

max(Ai  „,  A2„)  <  8s2.  (5.18) 

Boundary  conditions  (5. 16)  for  equation  (5. 15)  are  overdetermined  ones.  On  the  other  hand,  the 
QRM  is  well  suitable  for  solutions  of  such  problems.  Hence,  we  use  the  QRM  to  approximate 
functions  qn  k,  see  subsection  5.3  and  remark  5.1. 


5.2.  The  iterative  process 


Since  equations  (5. 15)  are  generated  by  the  Volterra-like  integral  differential  equation  (5.8),  it 
is  natural  to  solve  them  sequentially  starting  from  n  —  1 .  Let  £  e  (0,  1 )  be  a  sufficiently  small 
number.  Consider  a  function  /  (x)  e  C2  (R)  such  that 


xW  = 


l, 

between  0  and  1 

0, 


xe 

for x  e  (0,£)  U  (1  -f,  1), 
x  e  R\  (0,  1). 


(5.19) 


The  existence  of  such  functions  is  well  known  from  the  real  analysis  course.  We  choose  the 
first  guess  for  the  tail  function  V(j  (x)  as  in  subsection  6.2.  For  each  n  e  [1,1V],  we  perform  m 
iterations  with  respect  to  tails.  Hence,  we  obtain  three  finite  sequences  of  functions: 


{q„,k(x)} 


C N,m ) 

(»,k)=(l,l)> 


{V„,k  (x)} 


{ N,m ) 

(n, *)=(!, 1)> 


(n,k) 

r 


(x)} 


( N,m ) 

(«,*)=(!,!)’ 


xe  [0,  1], 


(5.20) 


Step  n(l\n  e  [1,1V].  Suppose  that  functions  qj(x)  e  H4  ( 0,  1)  and  V„_i(x)  e  C3[0,  1]  are 
constructed.  We  set  V,h\{x)  14- 1  (x).  In  particular,  V\t  \  (x)  :=  Vo(x).  Next,  using  the  QRM 

(subsection  5.3),  we  approximately  solve  equation  (5.15)  for  k  =  1  with  overdetermined 
boundary  conditions  (5.16)  and  find  the  function  qn  i  e  H4  (0,  1)  this  way  (remark  5.1). 
Hence,  by  the  embedding  theorem  qn  \  e  C3[0,  1],  Next,  we  find  the  approximation  £ *"’ 1 )  for 
the  unknown  coefficient  er  (x)  via  the  following  two  formulas: 


n- 1 

r„,  l  (x)  =  —hq,h  i  -  h  ^  qj  +  V„,  i ,  x  €  [0,  1], 

3=0 

e£n,1)(x)  =  1  +  r"  |  (x)  +  s2  [r’n ,j(x)]2  -  2snr’n  l(x),  x  e  [0,  1], 


(5.21) 

(5.22) 


pin,  1) 


(x) 


^"^(x),  if  £<'u,(x)  e  [1,  d],xe  [0,  1], 
1,  ifs^l)(x)  <  l,xe  [0, 1], 

d,  if  £,(n  l)(x)  >4xe[0,l]. 


(5.23) 
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Hence,  e  C1  [0,  1]  and  e^l)  e  C“  [0,  1],  Vcr  e  (0,  1).  Formulas  (5.21)  and  (5.22)  are 
obvious  discrete  analogues  of  formulas  (5. 1)  and  (5.6),  respectively.  Consider  the  function 

=  (1  -xW)  +  /(x)e',,'1)(x),  tel.  (5.24) 

It  can  be  easily  proven  that  e'/"' 1  *  (x)  e  [1,  d\  for  tel  This  implies  the  ellipticity  of  the 
operator  (2. 8)  for 

sr(x)  :=  er(,U)(x)  e  C  (K).  (5.25) 

Hence,  we  solve  the  forward  problem  (2.8),  (2.9)  and  (5.25)  for  s  :=  s.  Let  (x,  s)  be  the 
solution  of  this  problem  considered  in  theorem  3.1.  Then,  using  (5.7),  we  set  for  the  next  tail 
.  In  wnA  (x,  s)  -  In  iv o  (x,  s) 

Vn.2(x)  := - - -  (5.26) 

s 

Step  n{k\  n  e  [1,  TV] ,  k  e  [2,  m].  Suppose  that  functions  qj(x)  e  H4  ( 0,1)  and 
j  e  [0.  n  —  1],  Vnj t(x)  e  C:  1 " [0.  1]  are  constructed.  Using  the  QRM,  we  approximately 
solve  the  problem  (5.15)  and  (5.16).  This  gives  us  the  function  q„k  e  H4  (0,  1)  C  C3[0,  1], 
Next,  we  find  the  approximation  sl"'k)  e  C“[0,  1]  for  the  function  er  (x),  using  formulas  (5.21)- 
(5.23),  where  the  index  ‘U  is  replaced  with  the  index  ‘k’.  Suppose  now  that  n  <  N.  Then,  we 
construct  the  function  s'r(n'k)(x)  via  (5.24),  where  («,  1)  is  replaced  with  (n,  k)  and  solve  the 
forward  problem  (2.8)  and  (2.9)  with  er(x)  e4n  k) (x)  and  s  s.lf  k  <  in.  then  calculate 
the  new  tail  V„:k+i(x)  by  formula  (5.26)  where  indices  (n,  1)  and  (n,  2)  are  replaced  with 
(n,  k)  and  (n,  k+  1),  respectively.  If,  however,  n  <  N  and  k  —  m.  then  we  set 

V„(x)  :=  Vnm+i(x)  e  C3[0,  1],  q„(x)  :=  qnjn(x)  (5.27) 

and  go  to  the  step  (n  +  1)(1).  Let  n  —  N  and  k  =  m.  Hence,  sequences  (5.20)  are  constructed. 
Then,  we  set  er,gi0b(x)  :=  e  C“[0,  1]  and  stop  iterations.  The  function  Er.giobW  is  our 

final  solution  as  defined  in  definition  4.1. 

5.3.  The  quasi-reversibility  method  (QRM) 

The  QRM  was  first  proposed  in  [33].  We  refer  to,  e.g.,  [15,  18,  29]  and  references  cited 
there  for  some  further  developments.  The  QRM  can  solve  linear  ill-posed  boundary  value 
problems  for  many  PDEs,  including  problems  with  overdetermined  boundary  conditions.  It 
finds  least-squares  solutions.  Let 

n—  1 

t^n.kix)  — Ai  nh  ^2  qj  ^1  ,nVn.k  —'i 2 . . 

j= 0 

and  H„  k(x)  be  the  right-hand  side  of  equation  (5.15).  Then  anj,  e  C2[ 0,  1],  l/„j:  e  L?  (0,  1). 
Equation  (5. 15)  can  be  rewritten  as 

l -n.k  {^qn.kj  ■ —  t}n,k  =  Hn,k- 

Let  y  e  (0,  1)  be  the  regularization  parameter.  In  the  QRM,  we  minimize  the  following 
Tikhonov  functional  with  respect  to  the  function  qn  i t,  subject  to  the  boundary  conditions  (5. 16), 

Jy  (qn.ki  Hn.k)  —  II Ln,kQln,k)  "I-  T  ll^«.^llir4(0,l)'  (5.28) 

Remark  5.1.  There  is  a  peculiarity  of  this  technique  in  the  ID  case.  We  use  the  approximate 
least-squares  QRM  solutions  of  problem  (5.15)  and  (5.16).  On  the  other  hand,  it  seems  to 
be  at  the  first  glance  that  for  a  given  function  Vny  one  can  find  the  function  qnj:  via  solving 
equation  (5. 15)  with  only  two  out  of  three  boundary  conditions  in  (5. 16).  However,  our  attempt 
to  do  so  for  computationally  simulated  data  has  failed  to  produce  good  quality  results;  see 
remark  3.1  on  page  130  of  [30].  This  is  because  of  the  approximate  nature  of  our  mathematical 
model. 
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Lemma  5.1  (Carleman  estimate,  see  [29,  p  168]  and  [30]).  For  any  function  u  e  H 1  (0,  1) 

with  u  (0)  =  u'  (0)  =  0  and  for  any  parameter  k  f  \,  the  following  Carleman  estimate  holds: 

[  (u")2  e~2Xx  dx^—f  [8(u")2 +  \{u')2 +  k\u)2]e~2x-xdx.  (5.29) 

Jo  16  Jo 

Lemma  5.2  [30].  For  any  number  y  e  (0,  1)  and  for  each  function  F  e  Lid).  1),  there  exists 
a  unique  minimizer  uy  e  //4  (0,  1)  with  uy  (0)  =  u'y  (0)  =  0  of  the  functional  Jy  (u,  F) 
in  (5.28)  and 

\\ur  IIh4(o,i)  ^  Cy  FIIl2(o,d  > 

where  the  constant  C  >  0  is  independent  of  a,uk,  F,  uY,  y.  Let  ||  fln,t  || C|-0  i]  ^  ao ,  where 
fl0  =  const.  >  0,  Let  the  function  u  e  H 4  (0,  1),  Ln  ^  (u)  =  H  and  u  (0)  =  u'  (0)  =  0.  Then, 
there  exists  a  constant  K  =  K  (ao)  >  0  independent  ofF,  H,  uy,  u  and  ysuch  that 

\\uy  —  M|lif2(o,n  ^  K  (ll-F  ~  -^1Il2(0,1)  +  IImIIh4(0,i))  • 

The  proof  of  lemma  5.2  is  based  on  lemma  5.1,  the  Riesz  theorem  and  the  variational 
principle.  The  main  difference  between  lemma  5.1  and  Carleman  estimates  for  elliptic 
operators  in  the  n— D  (n  f  2)  case  is  that  we  now  have  the  integration  on  the  right-hand 
side  of  (5.29)  over  the  entire  interval  (0.  1)  rather  than  over  a  subdomain  of  the  original 
domain  in  n— D  [34].  The  absence  of  local  minima  of  the  functional  (5.28)  follows  from 
lemma  5.2.  It  can  be  derived  from  lemma  5.1  that  it  is  not  necessary  to  use  the  regularization 
term  y  \\q„,k  ||^4(0  1(  in  (5.28).  Still,  we  use  this  term  to  ensure  that  functions  e(rn'k)  eC'[0,  1], 
Although  lemma  5.1  implies  that  it  is  not  necessary  to  use  the  condition  q'nk(  1)  =  0,  we  still 
use  it  in  our  computations  to  improve  the  stability. 

5.4.  A  brief  scheme  of  the  numerical  method 

For  the  convenience  of  the  reader,  we  present  in  this  subsection  a  brief  summary  of  the  above 
numerical  method  by  outlining  its  main  steps. 

Step  1.  Calculate  the  Laplace  transform  (2.7)  Cg  =  (pis)  of  the  data  g  ( 1 )  in  (2.6). 

Step  2.  Using  the  function  (pis)  as  well  as  formulas  (5.3)  and  (5.5),  calculate  the  boundary 
data  i/ro0)  :=  q  (0,  s) ,  i/fi  (5)  :=  qx(0,s),s  e  [s,  s]  in  (5.9)  for  the  function 
q(x,s)  =  ds  {s~2  ln  [(w/wo)  (x,  i)]}.  Set  qo(x)  =  0.  Let  {q„(x)}^=1  be  functions 
generated  by  the  function  q(x,s),  as  defined  in  (5.12).  Use  (5.17)  to  calculate  the 
boundary  data  for  these  functions  at  [x  =  0} ,  fo,n  '■=  <7n(0),  '■=  4n  (0)- 

Step  3.  Calculate  the  first  tail  function  Vq(x)  as  in  subsection  6.2.  Set  Vi,i(a)  :=  Vq(x). 

Step  4.  Given  n  e  [l,N],fc  e  [1  ,m],  calculate  the  function  qn.k(x ).  To  do  this,  find  an 
approximate  solution  of  equation  (5.15)  with  boundary  conditions  (5.16),  using 
the  QRM,  via  the  minimization  of  the  functional  (5.28),  subject  to  boundary 
conditions  (5. 16).  As  soon  as  the  function  q,uk(x)  is  approximately  calculated,  find  the 
approximation  e,("’4)(x)  for  the  unknown  coefficient  er(x)  via  formulas  (5.21)-(5.23). 
Next,  find  the  function  e'fn,k)(x)  via  (5.24).  In  formulas  (5.21)-(5.24),  the  index  ‘1’ 
should  be  replaced  with  the  index  ‘k’ . 

Step  5.  Solve  the  forward  problem  (2.8)  and  (2.9),  with  s  :=  s  and  er(x)  :=  firn,k) (x).  Next, 
calculate  the  new  tail  function  Vnji+  \  (x)  using  (5.26),  where  indices  (n,  1)  and  (n,  2) 
should  be  replaced  respectively  with  (n,  k)  and  (n,  k  +  1 ). 

Step  6.  If  k  <  m,  then  go  to  step  4  with  k  :=  k  +  1.  If,  however,  k  =  m  and  n  <  N ,  then 

use  (5.27).  Next,  set  V„+ 1  1  (x)  :=  V„  (x)  and  go  to  step  4  with  n  :=  n  +  1,  k  :=  1. 

Finally,  if  k  —  m  and  n  —  N,  then  set  the  approximate  solution  of  the  above  CIP  as 
G.giob (x)  :=  £<rN'm'1  (x)  (see  definition  4.1)  and  stop  the  iterative  process. 
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6.  Convergence  analysis 


6.1.  The  exact  solution 


Using  again  the  concept  of  Tikhonov  for  ill-posed  problems  [6,  22,  25,  46],  we  assume  that 
there  exists  a  unique  exact  solution  e*(x)  of  our  CIP  with  the  ‘ideal’  noiseless  data  g*  ( t ) 
in  (2.6),  and  the  function  e*(x)  satisfies  conditions  (2.1)  and  (2.2).  Then  e*(x)  generates 
functions,  which  are  defined  similarly  with  the  ones  above,  q*0  =  0, 

r*(x,  s),  V*  (x,  s)  :=  r*  (x,  s) ,  q*{x,  s )  =  9sr*(x,  s),  q*(x),  ^0%,  «  e  [1,  N] .  (6.1) 


Similarly  with  (5. 15)  and  (5. 16),  we  obtain  from  (5.8)  and  (5.9), 


n—  1 


/  «— 1 


n— 1 


dxq*n  -  M.nh  J2  9 x4j  -  M.ndxV*  -  2A2,„  dxq*n  =  - A2,nh 2  ^  dxq*  +2 dx 


j= o 


\j= 0 


j= 0 


n—  1 


+  2 A2,ndxV*  |^/7  ^  dxq*J  -  A2,n  {dxV*Y  +  2 A2,ndxV* 
+  F|  „  (x,  /a,  h,  s) ,  xe  (0, 1), 


(6.2) 


q*n  (0)  =  .  9^:  (0)  =  rUl ,  9,<*  ( 1 )  =  0.  (6.3) 

In  (6.2),  Fn(x,n,h)  is  the  error  function  which  is  generated  by  averaging  (5.17)  of  the 
boundary  functions  t/tq  and  i/r*.  Also,  the  nonlinear  term  (/r,  h)  (dxq*)  is  included  in 
By  embedding  theorem  there  exists  a  constant  C2  Js  1  such  that 

ll/llc[o,i]  <  C2  11/llflUo.t) ,  V/  e  H1  (0, 1).  (6.4) 


By  the  above-mentioned  concept  of  Tikhonov  for  ill-posed  problems,  we  can  assume  that  we 
know  a  number  C*  >  0  such  that 


max 

ne[l,JV] 


\Wn 


,  C* 

Icqo.i]  ^ 


max 

H£[1,V] 


9, 


«  ll//4(0,l) 


<  C* 


(6.5) 


6.2.  Approximate  mathematical  model 

Following  definition  4.1,  we  introduce  in  this  subsection  an  approximate  mathematical  model. 
Assumptions  of  this  model  are  based  on  the  asymptotic  formulas  (5.10)  and  (5.11).  These 
assumptions  actually  mean  that  we  consider  only  the  first  term  of  the  asymptotic  behavior  of 
the  tail  function  V*  (x,  s)  when  the  pseudofrequency  s  — >  oo  and  truncate  the  rest.  We  call  s 
the  truncation  pseudofrequency.  A  similar  assumption  is  used  in,  e.g.  geometrical  optics.  Our 
approximate  mathematical  model  consists  of  the  following  two  assumptions. 

1 .  Let  the  parameter  s  >  1.  Then,  there  exists  a  function  a*  e  H4  (0,  1 )  such  that  the 
function  V *  (x,  ,v)  has  the  form 

V*(x,  s)  =  s~1a*(x),  Vs  >  s,  Vx  e  [0,  1],  (6.6) 

2.  The  following  equality  holds: 

s~1a*(x)  =  s_2[ln  w*(x,  s)  —  In  wo(x,  s)].  (6.7) 

Compare  (6.6)  and  (6.7)  with  (5.7)  and  (5.11).  Using  the  third  condition  (6.1)  and  (6.7),  we 
obtain 

q*  (x,  s)  =  —  s~2a*(x).  (6.8) 
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Consider  equation  (5.8)  for  the  pair  (q*(x,  s),V*  (x,s)) ,  s  e  [s,  s].  Then,  the  boundary 
conditions  (5.9)  are  valid  with  functions  t/tq  (s)  and  i j/*  ( s ).  Hence,  substituting  in  (5.8)  and  (5.9) 
s  s  and  using  (6.8),  we  find  that  the  function  a*  (x)  is  the  solution  of  the  following 
overdetermined  boundary  value  problem: 

d2a*(x)  =  0,  .16(0,1),  (6.9) 

a*  (0)  =  -;yV0*  (5),  dxa*  (0)  =  -s2f*  ( s ) ,  3 xa*  (1)  =  0.  (6. 10) 

In  (6. 10),  exact  boundary  data  at  x  —  0  are  used.  In  practice,  however,  we  have  the  non¬ 
exact  data  boundary  \l/o(s),  i// * (s ) ■  Hence,  consider  the  following  boundary  value  problem  for 
the  function  a(x): 

d2a(x)  =  0,x  e  (0,  1),  (6.11) 

a  (0)  =  -is Vo  (s) ,  3xa  (0)  =  -s2f\  (s) ,  dxa  (1)  =  0.  (6. 12) 

We  solve  problem  (6.11)  and  (6.12)  via  the  QRM  (see  remark  5.1).  Let  the  function  aY  e 
H 4  (0,  1)  be  the  unique  minimizer  of  the  following  analogue  of  the  QRM  functional  (5.28): 

JY(a)  =  ||a  lli2(o,n  +  Y  llallff4(o,i) ,  (6.13) 

subject  to  boundary  conditions  (6.12).  Because  of  (6.6)-(6. 13),  we  define  the  first  guess  for 
the  tail  function  as 

Vo  (x)  =  s~1ay(x).  (6. 1 4) 

Let  cr  e  (0,  1 )  be  the  level  of  the  error  in  the  boundary  data.  We  assume  that 

I  IMS)  -  ^o*(S)l  +  ItfrjCO  -  ^r*(s)  I  <  a.  (6.15) 

Remark  6.1.  Let  aY(x)  be  the  approximate  solution  of  problem  (6. 1 1) — (6. 13).  Then, 
substituting  (6.14)  into  (5.1)  and  (5.6)  at  s  =  is,  one  can  find  a  good  approximation  for 
the  exact  coefficient  s*(x).  Furthermore,  theorem  6.1  implies  that  all  functions  g(r"-k>  are  good 
approximations  for  e*  (x),  as  long  as  the  total  number  of  iterations  is  not  too  large.  This 
corresponds  well  with  (4.1).  Since  we  find  aY(x)  only  using  the  boundary  data,  this  means 
that  our  approximate  mathematical  model  is  indeed  a  good  one.  Hence,  we  can  stop  iterations 
on  any  function  s,n,k)  for  those  indices  («.  k),  which  are  ‘allowed’  by  theorem  6.1.  Next,  one 
can  use  the  adaptivity  procedure  to  refine  the  solution  (the  end  of  subsection  4.1).  This  was 
confirmed  numerically  in  tests  2  and  3  of  [7]  as  well  as  in  tests  2  and  3  in  section  4.16.2  of 
[6].  However,  if  not  using  the  adaptivity  for  refinement,  then,  quite  naturally,  one  needs  to  find 
an  optimal  iteration  number  to  stop;  see  figures  7.3,  7.5,  7.6  and  7.8  of  [5],  figures  3  and  4  of 
[10],  figures  6  of  [26]  and  figure  1(b)  of  [30]  (this  again  corresponds  well  with  definition  4.1). 
These  figures  can  also  be  found  in  chapters  3-5  of  [6],  along  with  objective  stopping  criteria 
for  iterations. 

Remark  6.2.  Because  of  the  approximate  nature  of  equalities  (6.6)  and  (6.7),  equation  (6.9) 
does  not  match  the  asymptotic  behavior  (5.10),  which  is  the  single  self-contradiction  of  this 
approximate  mathematical  model.  The  same  can  be  stated  about  all  other  versions  of  this 
method  in  the  above  cited  publications.  Nevertheless,  it  was  consistently  demonstrated  that 
this  numerical  method  works  well  for  both  computationally  simulated  and  experimental  data. 
Based  on  our  numerical  experience,  we  believe  that  this  is  because  of  two  factors:  (1)  the 
truncation  of  the  asymptotic  series  with  respect  to  \/s  at  s  -»  oo  in  (6.6)  is  reasonable  and  (2) 
the  procedure  of  updating  tail  functions. 
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Lemma  6.1.  For  each  value  of  the  regularization  parameter,  y  e  (0,  1 ),  there  exists  a  unique 
minimizer  ay  €  H 4  (0,  1)  of  the  functional  (6.13)  satisfying  the  boundary  conditions  (6.12). 
The  following  estimates  hold: 

IK  11^(0, D  <  Cs2y-l/2  (l^o  (5)1  +  IVn  (5)1), 

where  the  constant  C  >  0  is  independent  of  ay,  y ,  s,  i/ro(5),  (5).  Let  assumptions  1  and  2 

hold  and  the  function  a *  6  H 4  (0,  1)  be  the  solution  of  problems  (6.9)  and  (6.10).  Let  (6.15) 
hold.  Then,  there  exists  a  constant  K  >  0  independent  of  a*  (x),  ay(x),  y,  s,  fo  (5)  ,  \jf\  (s) 
such  that 

II Vo  -  VII^d  <  Ks(o  +  7/110*11^(0,!)). 

This  lemma  easily  follows  from  lemma  5.2,  (5.5)  and  (6.6)-(6. 14).  Uniqueness  within  the 
framework  of  this  approximate  mathematical  model  can  be  easily  derived;  see  lemma  6.6.2 
in  [6]  for  a  similar  result.  We  assume  below  that  the  above  exact  solution  e*  (x)  is  exactly  the 
same  as  the  one  within  the  framework  of  this  model,  i.e. 

e*(x)  =  r*xx  +  s2(r*x)2  -2sr*+  1,  se[5  s],  x  e  [0,  1];  (6.16) 

see  (5.1).  Here,  the  function  r*  is  defined  via  (5.6)  and  (5.7),  where  the  functions  q  and  V 
are  replaced  with  the  functions  q*  (x,  s)  and  V*  (x,  s ),  and  (6.6)— (6. 8)  hold.  Hence,  (6.2)-(6.4) 
and  (6. 16)  imply  the  following  analogues  of  the  discrete  formulas  (5.21)  and  (5.22): 

72—1 

r*  (x)  =  -hq* (x)  -hJ2  9*  M  +  U*  (x,  5)  +  F2,„  (x,  h.s),  x  e  [0,  1  ] ,  (6. 1 7) 

7=0 

s*(x)  =  1  +  d2r*(x)  +  sl  [dxr*(x)f  -  2 sndxr*(x)  +  FXn  (x,  h,  5) ,  x  e  [0,  1],  (6.18) 

One  can  prove  that 

3 

J2  Il^>lli2(0,i)  ^  C\s2(h  +  p-1),  (6.19) 

7=1 

where  the  constant  6)  >  0  is  independent  of  h,  p,  s.  We  assume  that 

Fj,n  =  0,j=  1,  2,  3;  tAo,„  =  if*n,  =  f*]n,  n  e  [1,  TV] .  (6.20) 

Therefore  by  (6.15)  the  error  is  ‘allowed’  only  in  numbers  (.v)  and  i//*  (,s).  We  point  out 
that  an  analogue  of  theorem  6.1  can  be  proved  very  similarly  for  the  realistic  case  (6. 19)  and 
when  |  foj,  —  i/r^  n  |  ,  |  i/xj  —  i jr*  n  \  <  a.  So,  condition  (6. 20)  is  introduced  for  brevity  only.  An 
additional  simplification  of  the  presentation  comes  from  balancing  the  error  parameters  h  and 
a  with  the  regularization  parameter  y  and  parameters  £  and  p  in  (5. 13)  and  (5. 19)  as 

a  =  7)7  =  =  p~l  h.  (6.21) 


6.3.  Approximate  global  convergence  theorem 

The  main  new  element  of  theorem  6. 1  is  that  we  iteratively  estimate  tails,  which  was  not  done 
in  [30].  To  do  this,  we  use  (3.25)  and  (3.26).  The  proof  of  theorem  6.1  has  some  similarities 
with  the  proof  of  theorem  6.7  of  [6]  for  the  3D  case.  However,  estimates  of  the  H2- norm 
of  the  QRM  solution  qnj:  in  [6]  are  obtained  in  a  subdomain  of  the  domain  fl.  In  turn,  this 
leads  to  the  estimate  of  the  accuracy  of  the  computed  target  coefficient  in  another  subdomain 
of  !T2  rather  than  in  the  entire  Q.  On  the  other  hand,  because  of  the  ID  case,  (5.29)  implies 
a  stronger  estimate  of  that  accuracy  in  the  whole  interval  (0,  1).  Still,  because  of  the  above 
similarity  with  the  proof  of  theorem  6.7  of  [6],  and  also  because  the  proof  of  theorem  6.1  is 
rather  technical,  we  only  outline  it  here. 
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Theorem  6.1  (Approximate  global  convergence).  Let  the  function  s*(x )  satisfying 
conditions  (2.1)  and  (2.2)  be  the  exact  solution  of  our  CIP  for  the  noiseless  data  g*  ( t ) 
in  (2.6).  Fix  the  truncation  pseudo  frequency  s  >  maxOT(r/,  |xo|),  1),  where  the  number 
T(d,  |xo|)  was  defined  in  theorem  3.1.  Let  two  assumptions  of  subsection  6.2  hold  and 
lla*ll//4(0 ,1)  ^  C*,  where  C*  is  the  constant  defined  in  (6.5).  Let  the  first  tail  function  Vq  (x) 
be  defined  via  (6.11)-(6.14).  In  addition,  let  conditions  (6.15),  (6.20)  and  (6.21)  hold.  Let 
B  =  B  ( xq ,  d,s)  >  1,  Ct,  C*  and  K  be  the  numbers  defined  in  theorem  3.2,  (6.4),  (6.5)  and 
lemma  6.1.  Define  B\  =  B\  (xq ,  d,  s)  =  max  (C2,  B,  C* ,  Ks).  Consider  the  iterative  process 
of  subsection  5.2.  Let  the  number  N  ^  2  be  independent  of  h.  Then  there  exists  a  constant 
D  —  D  (jco,  d,  s)  ^  32s2B\  such  that  if  for  a  constant  b  >  1  independent  of  x 0,  d,  s, 

h  e  (0,  h0),  h0=  bJNm+2,  (6.22) 

then  the  following  estimates  are  valid: 

jj(n-l)m+k+2fo 


ll^llcMo,!]  <  2C* 


(6.23) 

(6.24) 


Jn.k) 


r  11^2(0,1) 


sJ  D 


(n—  1  )m+k+ 1 


h, 


(6.25) 


\K,k  ~  9,V*||l2(o,i)  +  \\Vnk  ~  9;V*||l2(o,i)  <  D 


\(n—  l)m+k+l 


h, 


I  ^n,k  II  C[0,1]  ^  2Bu 


(6.26) 

(6.27) 


5  \  jfi  5^!  4  i  A\,nVU'k  2,n 

7=0 


^  18s2B]  :=  ag- 


(6.28) 


C[0,1] 


Here,  qn^  are  the  QRM  solutions  of  the  boundaiy  value  problems  (5.15)  and  (5.16).  In 
particular,  let  the  number  co  e  (0,  1)  be  defined  as 

In  b 


co  ■ 


In  b  +  (2Nm  +  2)  In  D 
Then  (6.25)  implies  the  Holder-like  estimate 

<  IC  :=  0, 


We^  -  e*\ 


£2(0,1) 


which  guarantees  the  approximate  global  convergence  property  of  the  level  6  of  this  iterative 
process  within  the  framework  of  the  approximate  mathematical  model  of  subsection  6.2 
(definition  4.1). 


Remark  6.3.  Since  the  number  N  of  5  subintervals  (,v,,  .v,_  1 )  is  assumed  to  be  independent 
of  the  partition  step  size  h,  then  theorem  6. 1  requires  the  length  of  the  total  .v  interval  [j,  5] 
covered  in  the  iterative  process  of  subsection  5.2  to  decrease  with  the  decrease  of  h.  This  seems 
to  be  a  natural  requirement.  Indeed,  if  the  number  s  —  y  would  be  independent  of  h,  then  this 
would  mean  the  increase  of  the  number  of  iteration  Nm  —  m  ■  (s  —  s)  //;.  On  the  other  hand,  the 
error  increases  with  the  iteration  number,  especially  for  nonlinear  ill-posed  problems.  Hence, 
the  number  of  iterations  Nm  is  one  of  the  regularization  parameters  here.  It  was  pointed  out  on 
pages  1 56—7  of  [22]  that  the  total  number  of  iterations  can  often  be  regarded  as  a  regularization 
parameter  in  the  theory  of  ill-posed  problems.  Two  other  regularization  parameters  are  s  and 
y.  Thus,  we  have  a  vector  (Nm,  k.  y  )  of  regularization  parameters. 
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Remark  6.4.  Another  argument  justifying  the  smallness  of  the  interval  [5,  s]  is  that  the  original 
equation  (5.8)  contains  Volterra-like  integrals  in  nonlinear  terms.  It  is  well  known  from  the 
standard  ODE  course  that  the  existence  of  the  solutions  of  Volterra-like  integral  equations  is 
guaranteed  only  on  sufficiently  small  intervals. 


Outline  of  the  proof  of  theorem  6.1.  To  avoid  using  new  notation,  below  qn,k(x)  means  the 
QRM  solution  of  problem  (5.15)  and  (5.16).  The  only  exception  is  when  we  subtract  (6.2) 
and  (6.3)  from  (5. 15)  and  (5. 16),  respectively.  Denote 

qn,k  =  qn,k  ~  q*n ,  ^  -  e*r,  Vn,k  =  V„k  -  V*. 

Consider  the  case  ( n ,  k)  —  (1,  1).  First,  by  lemma  6.1 


KiLto.i)  + 


v 


'U  lli2(0,l)  ^ 
Hence,  (3.25),  (6.4),  (6.22)  and  (6.29)  imply  that 


^  Bill  <  Dzh. 


Kil 


=  II V?  1  +  dxv* 


^  C2D  h  T  B\  2B\ 


(6.29) 


(6.30) 


llc[0,l]  II '1.1  1  x  11  C[0, 1] 

Estimates  (6.29)  and  (6.30)  establish  (6.26)  and  (6.27)  for  («,  k)  =  (1,  1).  Consider  now  the 
function  §11.  Subtracting  (6.2)  from  (5.15)  and  (6.3)  from  (5.16),  using  (6.20)  and  setting 
n  —  k  =  1 ,  we  obtain 

§U  +  (AuKi  +  2 A2a)  5j !  =  -  [Ahldxq\+A2,i  (V',  +  9vV*)]  Vu,  x  e  (0,  1) ,  (6.31) 

§1,1  (0)  =  §!,!  (0)  =  5ja  (1)  =  0.  (6.32) 

We  now  need  to  estimate  the  QRM  solution  q]  ,  of  problem  (6. 3 1 )  and  (6. 32).  By  (5. 1 8),  (6. 29) 
and  (6.30) 

| AhlVh  +  2A2,i  \  s?  8s2  (Bi  +  2)  <  \a.  (6.33) 

Next,  using  (3.25),  (5.18),  (6.5),  (6.29)  and  (6.30),  we  obtain 

\Ai,idxq\  +A2J(V;i  +  dxV*)\  <  8r(C*  +  3BQ  <  32 s2Bu 


\\[Ax,Aq\+A2,i{Vll  +  9^V*)]V1'il||  <  32 s2B\h  <  D~h. 


(6.34) 


We  can  take  D  ^  K  (ao),  where  numbers  a()  and  K  («o)  are  defined  in  (6.33)  and  lemma  5.2, 
respectively.  Hence,  lemma  5.2,  (6.31),  (6.32)  and  (6.34)  imply  that 

lkullfl*(0,i)  ^ 

which  proves  (6.23)  for  q\  \.  Next,  to  prove  (6.24),  we  use  (6.5),  (6.22)  and  (6.35), 
Ikullcqo.i]  ^  C2 ll^t.i  ll«2(o,t)  ^  C2lki,i ||rr2(o,i)  +  C2\\q\ ||h2(0,i)  <  2C* . 

Next,  subtract  (6. 18)  from  (5.22).  Note  that  by  (2.1)  and  (5.23) 

\e?'k)(x)\  =  \e(;'k\x)  -  <(x)|  <  | £<"■*>  (x)  -  <(x)|,  Vx  e  [0,  1], 

Hence,  using  (6.22),  (6.29),  (6.30),  (6.35)  and  (6.36),  we  obtain 

1  MHl2«>,1)  ^  (D3h2  +  Bit 0(1  +  f(3C*h  +  35 1  +  2)) 


(6.35) 


(6.36) 


^  Ts2(D3h2  +  Bih)  <  [4s2 Dh  ^  D2h , 


which  establishes  (6.25)  for  n  =  k  =  1.  Next,  using  (3.26),  (5.19),  (5.24)  and  (6.21),  we 
obtain 


V 


<  2 BD2h  +  ^2 Jd  <  D3h. 


„  „  V" 

!.2  II  i2C0, 1)  ^  II  1-2IIl2(0,1) 

Hence,  similarly  with  (6.30)  ||  V'j  2 1|  ^  2B\.  Thus,  (6.26)  and  (6.27)  are  established  for 
n  —  1 ,  k  —  2  The  rest  of  the  proof  can  be  done  similarly  using  mathematical  induction. 
We  need  (6.28)  to  estimate  norms  ||5/i,r  ||H2(0  1(  via  lemma  5.2.  The  estimate  (6.28)  is  proved 
using  (5. 18),  (6.22)  and  (6.24).  □ 
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Figure  1.  A  schematic  diagram  of  data  collection  by  the  Forward  Looking  Radar  of  the  US  Army 
Research  Laboratory. 


7.  Experimental  setup  and  data  pre-processing 


7.1.  Data  collection 


A  schematic  diagram  of  the  data  collection  by  the  Forward  Looking  Radar  is  depicted  in 
figure  1 .  Time-resolved  electric  pulses  are  emitted  by  two  sources  installed  on  the  radar.  Only 
one  component  of  the  electric  field  is  both  originated  and  measured  in  the  backscattering 
regime.  The  data  are  collected  by  16  detectors  with  a  step  size  in  time  of  0.133  ns.  Only 
shallow  targets  placed  both  below  and  above  the  ground  can  be  detected.  The  depth  of  the 
upper  part  of  a  shallow  underground  target  is  a  few  centimeters.  A  ground  positioning  system 
(GPS)  provides  the  distance  between  the  radar  and  a  point  on  the  ground  located  above  that 
target.  The  error  in  the  latter  is  of  a  few  centimeters.  Time-resolved  voltages  of  backreflected 
signals  are  integrated  over  radar/target  distances  between  20  and  8  m,  and  they  are  also 
averaged  with  respect  to  both  source  positions  and  with  respect  to  readings  of  16  detectors. 
Since  the  radar/target  distance  is  known,  it  is  approximately  known  which  part  of  the  measured 
time-resolved  signal  corresponds  to  the  reflections  from  that  target;  see  figure  1.  However, 
clutter  obviously  obscures  a  significant  part  of  that  signal.  For  any  target  of  interest,  only  a 
single  time-dependent  curve  can  be  extracted  from  the  vast  amount  of  data;  see  samples  in 
figures  4(b),  (d)  and  (f).  This  is  the  curve  we  have  worked  with  in  each  of  the  available  five 
cases  of  experimental  data. 

Since  the  radar/target  distance  is  provided  by  the  GPS  with  a  good  accuracy,  geometrical 
parameters  of  targets,  including  their  depths,  are  not  of  interest  here.  The  main  goal  of  our 
work  was  to  calculate  ratios  R  of  dielectric  constants 


Sy  (target) 
r  (bckgr)  ’ 


where  sr  (bckgr)  is  the  dielectric  constant  of  the  background  medium.  If  sr  (bckgr)  is  known, 
then  (7.1)  enables  us  to  calculate  er  (target).  If  a  target  is  located  above  the  ground,  then  er 
(bckgr)  =  sr  (air)  =  1.  Since  targets  can  be  mixtures  of  constituent  materials,  the  sr  (target)  is 
a  certain  weighted  average  of  dielectric  constants  of  these  materials.  We  image  the  ratio  (7. 1) 
rather  than  the  function  sr(x)  itself,  since  (2. 1)  requires  that  er  (bckgr)  should  have  a  constant 
value  outside  of  our  domain  of  interest  x  e  (0,  1).  The  latter  was  true  only  for  targets  located 
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above  the  ground.  However,  in  the  case  of  underground  targets,  this  condition  is  invalid  because 
of  the  air/ground  interface.  Below  ‘1’  stands  for  ‘1  m’. 

A  separate  question  is  about  the  meaning  of  the  dielectric  constant  in  metallic  targets. 
Comparison  of  figures  2(a)  and  (b)  shows  that  the  electric  held  reflected  from  the  metallic  target 
is  about  the  same  as  that  reflected  from  the  dielectric  target  with  a  large  value  of  the  dielectric 
constant.  Next,  we  simulated  the  data  for  four  large  inclusion/background  contrasts  by  solving 
(2.8)  and  (2.9).  Figure  2(c)  shows  the  computed  functions  w  (0,  s )  —  Wq  (0,  s),  s  e  [1.  12] 
for  the  case  of  a  single  inclusion  embedded  in  the  interval  (0,  1)  with  four  different  values  of 
the  function  sr  =  10,  20,  30,  40  in  this  inclusion;  see  figure  3(a)  for  the  geometry.  One  can 
observe  that  these  curves  do  not  change  much  with  a  change  of  sr  e  [10,  30].  Furthermore, 
curves  for  sr  —  30,  40  are  almost  the  same.  Therefore,  based  on  figures  2(a)-(c),  we  choose 
interpretation  (7.2)  of  the  dielectric  constant  in  metals.  The  physical  meaning  of  (7.2)  is 
unknown,  and  we  call  it  appearing  dielectric  constant  for  metals, 

£r (metal)  e  [10,  30],  (7.2) 


7.2.  Main  uncertainties  in  the  experimental  data 

To  figure  out  what  kind  of  ideal  data  one  should  expect  for  the  case  of  one  target  only,  we 
performed  computational  simulations  via  solving  the  forward  problem  (2.3)  and  (2.4)  for  the 
case  of  one  target  and  for  the  source  position 

*o  =  -l.  (7.3) 

In  data  simulation,  we  replaced  R  in  (2.3)  with  the  interval  x  e  (—6,  6)  and  have  considered 
zero  Dirichlet  boundary  conditions  for  rather  small  times, 

u  (—6,  t)  —  u  (6,  f)  =  0,  t  e  (0,  4).  (7.4) 

Condition  (7.4)  is  imposed  because  the  wave  front  originated  at  the  source  position  (7.3)  does 
not  reach  points  x  =  ±6  for  1  e  (0,4).  The  structure  of  the  medium  and  the  computed 
function  u  (0,  t)  :=  g  ( t )  are  depicted  in  figures  3(a)  and  (b),  respectively.  In  the  case 
sr  (target)  =  const,  e  (0,1)  the  function  g  (f )  looks  similar  (not  shown),  except  that  its 
peak  points  downwards.  Note  that  the  constant  background  in  figures  3(a)  and  (b)  corresponds 
to  the  fundamental  solution  of  the  ID  wave  equation  vtl  =  vxx, 

Uq  (X,  t )  =  \H  (t  —  \x  —  Xol), 

where  H  is  the  Heaviside  function.  Figures  4(b),  (d)  and  (f)  show  the  experimental  data  for 
different  targets.  A  visual  comparison  of  figures  4(b),  (d)  and  (f)  with  figure  3(b),  confirms  the 
above  mentioned  (section  1)  substantial  discrepancy  between  computationally  simulated  and 
experimental  data.  This  discrepancy  is  the  main  challenge  of  working  with  these  data. 

In  addition  to  the  above  misfit  and  the  3D  reality  versus  only  a  single  curve  for  each  target, 
there  were  some  other  uncertainties  here  as  well.  The  most  significant  uncertainties  were  as 
follows. 

(1)  The  reference  signal  was  not  measured. 

(2)  The  direction  of  the  incident  plane  wave  was  oblique  to  the  ground  rather  than  orthogonal; 
see  figures  4(a),  (c)  and  (e). 

(3)  Units  for  the  amplitude  of  experimental  data  were  unknown. 

(4)  The  location  of  the  point  source  xq  was  unknown.  Thus,  (7.3)  is  an  intuitive  choice. 

(5)  The  time  moment  t  =  0  on  the  data  was  unknown. 

(6)  The  background  was  heterogeneous  due  to  clutter. 
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Figure  2.  Appearing  dielectric  constant  for  metals  for  computationally  simulated  data  at  a  single 
frequency,  (a)  The  amplitude  of  the  electric  field  reflected  from  a  dielectric  target  with  sr  =  10  in  it. 
(b)  The  amplitude  of  the  electric  field  reflected  from  a  metallic  target  of  the  same  frequency  as  that  of 
(a).  Reflected  fields  (a)  and  (b)  are  very  similar,  (c)  Computed  functions  (w  —  wo)  (0,  s),  s  e  [1 ,  12] 
for  the  case  when  one  inclusion  is  embedded  in  the  interval  (0,  1)  with  different  values  of  the 
function  sr;  see  figure  3(a)  for  the  location  of  the  inclusion.  From  top  to  bottom  sr  =  10,  20,  30,  40. 
Curves  for  sr  =  30,  40  almost  coincide.  Hence,  the  function  w  (0,  5)  does  not  change  much  with 
an  increase  in  the  inclusion/background  contrast  from  10  to  40.  (a)-(c)  justify  definition  (7.2)  of 
the  appearing  dielectric  constant  for  metals. 
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Figure  3.  Computationally  simulated  data  u  (0,  t ),  where  u  (x,  t )  is  the  solution  of  problem  (2.1) 
and  (2.2)  with  the  source  location  Jto  =  —1.  Boundary  conditions  (7.4)  were  used,  since  the  wave 
front  did  not  yet  reach  points  x  =  d=6  for  times  t  e  (0,  4).  (a)  The  function  £r(x), x  €  (0,  1).  Also, 
er(jt)  =  1  for  x  $.  (0,  1).  (b)  The  function  u  (0,  t )  in  time  T  =  [0,  3.5]. 
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Remark  7.1.  It  is  because  of  items  4  and  5  that  it  was  impossible  to  image  locations  of  targets 
correctly. 

At  the  same  time,  the  UNCC/ChalmersGU  team  had  the  following  two  pieces  of  the 
information  in  advance. 

(1 )  The  presence  of  a  single  target  for  each  data  set. 

(2)  It  was  known  whether  the  target  was  located  above  or  below  the  ground. 


7.3.  Data  pre-processing 

We  need  to  pre-process  the  experimental  data  in  such  a  way  that  the  resulting  time- 
resolved  curves  would  look  similar  to  figure  3(b),  since  this  would  fit,  at  least  somehow, 
our  mathematical  model.  If  the  target  is  located  above  the  ground,  then  er  (target)  >  1 ,  since 
e,(bckgr)  =  er  (air)  =  1  in  this  case.  Figure  3(b)  indicates  that  one  should  select  on  the 
experimental  curve  only  one  downward-looking  peak  in  this  case.  However,  if  the  target  is 
buried  in  the  ground,  then  there  could  be  any  relation  between  e,(  target)  and  e,  (bckgr).  Hence, 
based  on  figures  3(a)  and  (b)  as  well  as  their  analogues  for  the  case  e,  (target)  <  e,  (bckgr)  (not 
shown),  we  selected  on  each  experimental  curve  the  earliest  peak  of  the  largest  amplitude  out 
of  all  other  peaks.  The  rest  peak  of  each  curve  was  set  to  zero.  More  precisely,  our  selection 
of  that  peak  was  as  follows:  this  should  be  the  earliest  peak  of  the  largest  amplitude 

out  of  Pea^s  f°r  a  target  buried  in  the  ground, 

all  downward-looking  peaks  for  a  target  above  the  ground. 

We  assigned  on  each  experimental  curve  the  time  zero  {t  =  0}  to  be  such  a  point  on  the  time 
axis,  which  is  1  ns  off  to  the  left  from  the  beginning  of  the  selected  peak.  Next,  we  multiplied 
the  resulting  data  by  the  scaling  factor  (below)  and  have  regarded  the  resulting  curve  as  the 
pre-processed  data.  Figure  5(a)  displays  the  pre-processed  data  for  the  case  of  figure  4(b). 
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Bush  (Clutter,  Test  2) 


Target  is  desert  bush  seen  along  the  road. 


(a) 


Wood  Stake  (Test  3) 


Target  is  wood  stake  on  the  ground. 

(C) 


(d) 


Figure  4.  Targets  and  experimental  data  for  tests  1-3.  The  ground  is  dry  sand  with  sr  e  [3,  5] 
[44].  Because  of  the  blind  study  case,  (a),  (c)  and  (e)  were  released  to  the  UNCC/ChalmersGU 
team  only  after  computations  were  made,  (a)  Bush  standing  along  the  road  (clutter),  (b)  Scaled 
experimental  data  for  (a)  (subsection  7.3).  The  horizontal  axis  is  time  in  nanoseconds  with  a  time 
step  of  0.133  ns.  The  vertical  axis  is  the  amplitude  of  the  measured  voltage,  (c)  Wood  stake,  (d) 
Scaled  experimental  data  for  (c).  (e)  Metal  box  buried  in  the  dry  sand,  (f)  Scaled  experimental 
data  for  (e).  The  amplitude  of  the  largest  downward-looking  peak  on  (f)  is  0.0072,  whereas  the 
amplitude  of  the  largest  upward-looking  peak  is  0.007.  Therefore,  the  downward-looking  peak 
was  used  in  our  data  pre-processing.  A  huge  misfit  between  experimental  and  computationally 
simulated  data  is  evident:  compare  highly  oscillatory  curves  of  (b),  (d)  and  (f)  with  figure  3(b). 
Waveforms  of  (b),  (d)  and  (f)  show  why  the  radar  detection  and  discrimination  problem  is  so 
challenging.  One  can  see  three  very  different  types  of  targets,  yet  their  signatures  are  very  similar. 


Figure  5(b)  shows  superimposed  pre-processed  curves  for  all  five  cases  of  experimental  data 
we  possess. 

The  amplitude  of  the  time -resolved  signal  for  each  case  was  of  the  order  of  105.  This  is  well 
above  the  amplitude  of  figure  3(b).  Thus,  all  pre-processed  data  were  multiplied  by  the  scaling 
number  SN  —  1  O'  7.  Here  we  show  how  we  have  chosen  the  number  SN  =  1 0  7.  Consider 
the  pre-processed  signal  for  the  bush  and  multiply  it  by  10  7 .  We  obtain  the  signal  depicted 
in  figure  5(a).  Next,  calculate  the  Laplace  transform  for  s  e  [1, 5]  of  two  signals:  (1)  that  of 
figure  5(a)  and  (2)  that  of  the  computationally  simulated  data  of  figure  3(b).  Superimposed 
graphs  of  the  function  w(x,  s)  =  w  (0,  s)  —  wq  (0,  s),  s  e  [1,5]  for  both  these  cases  are 
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Flush  Buried  Metal  Box  (Test  4) 

Radar  line  of  sight 


Air 


Dry  sand:  f*=3-5 


scmy=y/40c 

< - >  v 

40  cm 


Target  is  flush  buried  metal  box. 

(e) 

Figure  4.  (Continued.) 


(f) 


Figure  5.  The  horizontal  axis  is  time  in  nanoseconds  with  a  time  step  size  of  0.133  ns.  (a)  The 
pre-processed  experimental  data  for  the  bush  standing  along  the  road  (clutter);  see  figures  4(a) 
and  (b).  (b)  The  pre-processed  experimental  data  for  all  five  cases  (superimposed).  In  terms  of  our 
mathematical  model,  we  use  these  curves  as  functions  u  (0,  t )  —  uq  (0,  t ). 


displayed  in  figure  6.  One  can  observe  in  this  figure  that  minimal  and  maximal  values  of  both 
curves  are  approximately  the  same. 

In  fact,  we  initially  tried  three  scaling  numbers:  SN\  —  1 0  ,  SN2  =  1  (0  7  and  SW3  =  1 0  s. 

In  the  end  we  have  chosen  SN2  SN  —  10-7  out  of  these  three.  We  have  made  this  choice 
because  of  our  observation  that  this  was  the  only  case  out  of  three  when  minimal  and  maximal 
values  of  functions  w(x,  s),  s  e  [1,5]  of  both  above  curves:  that  of  the  Laplace  transform 
of  the  function  depicted  in  figure  3(b)  and  that  of  the  Laplace  transform  of  the  function  of 
figure  5(a),  were  approximately  the  same,  see  figure  6.  On  the  other  hand,  for  57V [  and  SN3  the 
minimal  and  maximal  values  were  quite  different  from  those  of  the  Laplace  transform  of  the 
function  of  figure  3(b).  As  soon  as  SN  =  I  0-7  was  chosen,  using  only  the  data  for  the  bush,  we 
multiplied  the  other  four  pre-processed  signals  by  10~7,  and  observed  similar  behavior  of  the 
minimal  and  maximal  values  of  the  Laplace  transforms  of  the  resulting  four  functions.  So,  the 
graphs  of  figure  5(b)  represent  five  available  pre-processed  signals  being  multiplied  by  1 0  7 
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Figure  6.  Graphs  of  the  function  w  (x,  s)  =  w  (0,  j)-«Jo(0,j),ie  [1,5]  for  the  Laplace  transform 
of  the  computationally  simulated  data  for  the  function  er(x )  displayed  in  figure  3(a)  and  for  the 
signal  of  the  bush  of  figure  5(a).  Recall  that  graphs  of  figures  5(a)  and  (b)  were  obtained  after 
multiplying  the  pre-processed  signals  by  the  scaling  number  SN  =  1 0  7  One  can  observe  that 
minimal  and  maximal  values  of  the  function  w(x,  s)  are  approximately  the  same  for  both  cases.  A 
similar  observation  was  in  place  for  four  other  cases  of  figure  5(b). 


In  the  case  of  the  upward-looking  peak  of  figure  5(b),  we  compared  the  function  \w(x,  s) 
for  it  with  the  function  w(x,  s )  for  the  simulated  data  above.  Graphs  of  figures  4(b),  (d)  and 
(f)  were  obtained  from  original  experimental  curves  by  multiplying  them  by  SN  —  10~7.  In 
addition,  we  conducted  a  limited  sensitivity  study  with  respect  to  the  scaling  factor;  see  test  2 
in  sections  8  and  9. 


7.4.  Functions  t/t0  (s)  and  i//|  (.S') 


Let  'git)  be  the  pre-processed  experimental  data  for  any  of  our  targets.  Based  on  figure  5,  we 
calculated  the  Laplace  transform  (2.7 )cp(s)  =  C  (g)  by  integration  over  the  interval  t  e  (0,  2). 
Next,  we  calculated  the  derivative  (p'(s)  as 


<p'(s)  =  ~  [ 

Jo 


g(t)te  *'df  +  9sw0(0,  s), 


(7.5) 


where  the  function  w o  (0,  ,v)  is  defined  in  (2.1 1).  It  was  observed  in  computational  simulations 
of  [30]  that  the  function  (p(s)  has  the  best  sensitivity  to  the  presence  of  inclusions  for 
s  e  [0.5,  1.2].  Nevertheless,  we  observed  that  the  larger  interval  s  e  [1,  12]  provides  better 
quality  images  for  simulated  data.  The  function  cp(s )  was  computed  accurately  for  the  entire 
interval  s  e  [1,  12]  with  the  step  size  As  =  0.05  in  the  .v-direction  However,  because  of  their 
dependence  on  the  derivative  tp\s)  in  (7.5),  functions  i I/q(s)  =  q( 0,  s),  i —  qx( 0,  j)  have 
oscillated  for  s  e  [3,  12].  On  the  other  hand,  our  testing  of  computationally  simulated  data  has 
shown  that  oscillations  should  not  be  present. 
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Hence,  we  pre-processed  the  function  fais)  as  follows.  First,  we  calculated  fais)  for 
s  e  [1,  2.5]  using  the  Laplace  transform  of  the  data,  as  in  (7.5),  with  the  0.5  step  size  with 
respect  to  Next,  we  set  i^o  (12)  :=  0.025  •  i/'i)  (2.5),  then  we  linearly  interpolated  the  plane 
(s,  i^o)  between  points  (2.5,  fo  (2.5))  and  (12,  i^o  (12))  and  similarly  for  fiis).  Having 
functions  fois)  and  i/qU)  for  our  pre-processed  experimental  data,  we  calculated  numbers 
i/fo,/i  and  i//]  „  in  (5. 17),  which  are  used  as  boundary  conditions  for  the  QRM. 


8.  Imaging  results 

A  detailed  description  of  the  numerical  implementation  of  the  algorithm  of  subsections  5.2 
and  5.3  can  be  found  in  [30],  We  now  briefly  outline  the  main  elements  as  well  as  some 
differences  with  [30].  Unlike  [30],  where  only  synthetic  data  were  used,  the  values  q„  (1)  were 
unknown  to  us  now  (naturally).  But  this  turned  out  not  to  be  a  problem.  Indeed,  our  numerical 
experiments  with  computationally  simulated  data  have  shown  to  us  that  knowledge  of  the 
numbers  qn  (1)  affects  only  the  accuracy  of  the  image  of  the  location  of  a  target.  However,  it 
does  not  affect  the  accuracy  of  the  reconstruction  of  the  target/background  contrast  in  sr,  which 
is  our  main  goal  here;  see  remark  7.1.  Thus,  unlike  [30],  we  did  not  assign  any  values  to  q„  ( 1 ) . 
The  initial  tail  function  was  taken  as  in  (6. 14).  Although  theorem  6.1  guarantees  approximate 
global  convergence  only  for  the  case  of  a  small  s  interval,  our  computational  experience  tells 
us  that  the  interval  s  e  [1,  12]  is  an  optimal  one.  We  attribute  this  to  the  well-known  fact  that 
convergence  estimates  are  routinely  much  more  pessimistic  than  computational  results.  This 
is  because  constants  in  convergence  theorems  usually  are  largely  overestimated.  Similarly, 
although  the  above  theory  works  only  for  the  case  sr(x)  f  1,  this  did  not  prevent  us  from 
computing  one  case  with  R  <  1  in  (7.1). 

We  regard  R  :=  R(x)  in  (7.1)  as  an  x-dependent  function.  With  respect  to  the  results  of 
this  section,  sr(x )  in  (2.8)  was  replaced  with  R(x).  Thus,  it  is  R(x)  which  was  computed  by  the 
above  algorithm.  Let  Rcomp  (x)  be  the  computed  coefficient  R(x).  Then  we  define  the  computed 
target/background  contrast  as  R  =  max  Rcmnv(x)  in  the  case  when  max  Rcmu]1(x)  >  1,  and  as 
R  =  minRcomp(x)  in  the  case  when  max  Raimp(x)  f  1.  We  set  sr  (target)  :=  R  ■  er  (bckgr). 

For  each  test,  we  also  computed  two  curves:  the  Laplace  transform  of  the  pre-processed 
experimental  data  and  the  function  w(0,  .v),  where  w  (0,  .v)  is  the  solution  of  problem  (2.8) 
and  (2.9)  with  e,  (x)  :=  /AompOO-  The  interval  .v  e  [1,  12]  with  the  step  size  h  —  0.5  was  used 
for  the  latter.  These  two  curves  were  very  close  to  each  other  for  all  tests. 

Test  1.  Computationally  simulated  data.  First,  we  verify  that  our  algorithm  provides  an 
accurate  target/background  contrast  for  computationally  simulated  data.  We  image  the  structure 
depicted  in  figure  3(a).  Figure  7(a)  displays  the  resulting  image.  Recall  that  we  are  interested  in 
accurate  imaging  of  target/background  contrasts  rather  than  in  accurate  imaging  of  locations 
of  targets  (subsection  7.1).  The  imaged  target/background  contrast  is  3.8,  whereas  the  real 
contrast  is  4.  Thus,  the  imaged  contrast  is  quite  accurate.  This  gives  us  a  hope  that  contrasts 
for  experimental  data  are  also  computed  accurately. 

Test  2.  The  image  of  a  bush  (see  figures  4(a)  and  (b),  5(a)  and  (b)).  This  was  the  most  difficult 
case  because  the  target  was  a  highly  heterogeneous  one.  Moreover,  the  maximum  of  the 
modulus  of  the  amplitude  of  the  pre-processed  signal  for  this  target  exceeds  these  values  for 
other  targets  by  a  factor  of  2.57.  Figures  7(b)  and  (c)  display  the  resulting  image  and  the  above 
superimposed  curves,  respectively.  Only  a  small  difference  between  the  curves  of  figure  7(c) 
is  observed. 

We  conducted  a  limited  sensitivity  study  with  respect  to  the  choice  of  the  scaling  number 
SN.  We  varied  SN  by  20%.  Figure  7(d)  displays  the  computed  function  R(x)  for  the  data  for  the 
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Figure  7.  Resulting  images,  (a)  The  computed  image  for  computationally  simulated  data  (test  1). 
Solid  and  dashed  lines  are  true  and  computed  images,  respectively.  Recall  that  we  are  interested 
in  accurate  imaging  of  target/background  contrasts  rather  than  in  accurate  imaging  of  locations  of 
targets  (subsection  7.1).  The  computed  target/background  contrast  is  3.8,  which  is  5%  error,  (b) 
The  computed  image  of  the  bush;  see  figure  4(a).  Here  R  =  6.4,  which  is  in  the  range  of  tabulated 
values  sr  e  [3,  20]  [20].  (c)  The  solid  line  is  the  Laplace  transform  of  the  function  depicted  in 
figure  5(a).  The  dashed  line  is  the  function  w(0,  s )  —  u>o(0, s),  where  w  (0,  s )  is  the  solution  of 
problem  (2.8)  and  (2.9)  with  £r(x)  :=  RComp(*)  of  (b).  (d)  Results  of  a  limited  sensitivity  study 
with  respect  to  the  choice  of  the  scaling  number  SN  =  0.8  x  10~7,  10-7,  1.2  x  10~7  for  the  case 
of  the  bush  (recall  that  we  chose  SN  =  10-7  for  all  our  tests).  One  can  observe  that  the  computed 
value  of  R  varies  as  R  =  4.5,  6.4,  8.6  All  these  three  values  are  within  tabulated  limits;  see  table  1. 
A  similar  observation  took  place  for  four  other  targets  we  worked  with,  (e)  The  computed  image  of 
the  wood  stake  (see  figure  4(c)).  Here  R  =  3.8,  which  is  in  the  range  of  tabulated  values  sr  G  [2,  6] 
(see  [44]).  (f)  The  computed  image  of  the  metal  box,  R  =  3.8. 


bush  obtained  for  SN  —  0.8  x  10”7,  10-7,  1.2  x  10~7.  Recall  that  we  chose  SN  =  10~7  for  all 
five  sets  of  experimental  data  we  possess  (subsection  7.3).  One  can  observe  that  corresponding 
values  of  the  computed  target/background  contrast  R  are  R  —  4.5,  6.4,  8.6.  The  value  R  —  6.4 
is  for  SN  —  1 0  7  and  is,  therefore,  the  same  as  that  of  figure  7(b).  All  these  three  values  are 
within  tabulated  limits,  see  table  1.  The  same  study  was  conducted  for  four  other  data  sets  we 
worked  with.  And  the  same  observation  was  in  place:  the  above  three  values  of  SN  resulted 
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(e)  (f) 


Figure  7.  (Continued.) 


Table  1.  Blindly  computed  values  R  for  five  targets,  £r, comp  (target)  :=  R  ■  sr  (bckgr).  Here  ‘a/b’ 
means  ‘above/below  the  ground’.  For  ‘a’  and  ‘b’  the  background  is  air  and  dry  sand,  respectively. 


Target 

a/b 

R 

e, (bckgr) 

£r,comp  (target) 

£r,pubi  (target) 

Figure  6(a) 

- 

3.8 

1 

3.8 

4,  figure  6(a) 

Bush 

a 

6.4 

1 

6.4 

G  [3,  20],  see  [20] 

Wood  stake 

a 

3.8 

1 

3.8 

G  [2,  6],  see  [44] 

Metal  box 

b 

3.8 

6  [3,  5],  see  [44] 

G  [11.4,  19] 

G  [10,  30],  see  (7.2) 

Metal  cylinder 

b 

4.3 

G  [3,  5],  see  [44] 

G  [12.9,21.4] 

G  [10,  30],  see  (7.2) 

Plastic  cylinder 

b 

0.28 

e  [3,  5],  see  [44] 

e  [0.84,  1.4] 

1.2,  see  [45] 

in  values  of  R ,  which  were  within  tabulated  limits  (corresponding  images  are  not  shown  for 
brevity);  also  see  section  9  for  a  discussion  of  this  issue. 

Test  3.  The  image  of  a  wood  stake  (see  figures  4(c)  and  5(b)).  The  computed  image  is  displayed 
in  figure  7(e). 

Test  4.  The  image  of  a  metal  box  (see  figures  4(e)  and  5(b)).  The  computed  image  is  displayed 
in  figure  7(f). 

Since  we  had  total  five  sets  of  data  in  our  possession,  we  also  imaged  two  more  cases: 
plastic  cylinder  and  metal  box,  both  buried  in  soil.  Dielectric  constants  were  not  measured 
when  the  data  were  collected.  Therefore,  we  compared  computed  values  of  dielectric  constants 
with  those  listed  in  tables  [44,  45].  Note  that  these  tables  often  provide  a  range  of  values  rather 
than  exact  numbers.  The  soil  was  dry  sand,  where  the  dielectric  constant  varies  between  3  and 
5  [44].  Denote  £,.pubi  (target)  the  published  value  of  the  dielectric  constant  for  non-metallic 
targets.  By  (7.2)  we  regard  £r,Pubi(target)  e  [10,  30]  for  metals.  The  case  of  bush  (vegetation) 
is  not  listed  in  [44,  45].  Hence,  we  took  er>pubi (target)  for  vegetation  from  figures  2  and  3  of 
[20].  Table  1  summarizes  our  results. 


9.  Discussion 

Since  dielectric  constants  were  not  measured  in  experiments,  the  maximum  that  can  be  done 
is  to  compare  computational  results  with  tabulated  values  (also,  see  the  last  paragraph  of  this 
section).  Therefore,  the  most  important  conclusion  from  table  1  is  that  computed  dielectric 
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constants  of  targets  fall  well  within  tabulated  limits  for  five  out  of  five  available  data  sets.  This 
was  achieved  regardless  of  the  significant  limitations  of  our  data.  Naturally,  these  limitations 
have  resulted  in  a  large  and  yet  unknown  noisy  component  in  the  data,  which  we  have  used  as 
inputs  for  our  algorithm.  Thus,  we  conclude  that  the  above  results  point  towards  the  robustness 
of  our  algorithm. 

Because  of  the  discussion  of  this  section,  it  is  worth  summarizing  now  the  severe 
limitations  of  our  data. 

(1)  Time-resolved  signals  were  integrated  and  averaged  over  radar/target  distances  between 
20  and  8  m  and  then  averaged  with  respect  to  two  sources  and  16  detectors. 

(2)  A  huge  misfit  between  experimentally  measured  and  computationally  simulated  data. 

(3)  The  3D  reality  versus  the  necessity  to  model  the  process  as  ID:  because  we  had  only  one 
time -resolved  curve  for  each  target. 

(4)  The  need  to  use  only  one  ID  hyperbolic  PDE  for  modeling  the  process  instead  of  the  full 
Maxwell  system:  because  only  one  component  of  the  electric  field  was  measured. 

(5)  A  number  of  other  uncertainties  listed  in  subsection  7.2. 

It  is  unlikely  that  an  accurate  image  of  a  3D  target  can  be  calculated  using  only  a  single 
time-dependent  curve  (item  3).  Therefore,  given  the  above  limitations,  we  had  only  a  very 
limited  goal  in  this  study.  More  precisely,  this  goal  was  to  calculate  a  single  number  for  each 
target,  which  would  characterize  the  target/background  ratio  of  dielectric  constants,  and,  at  the 
same  time,  would  fall  within  tabulated  limits.  That  number  is  R.  We  point  out  that  the  number 
R  is  still  computed  by  first  computing  the  function  R(x)  in  (7. 1). 

The  data  pre-processing  procedure,  which  was  similar  to  that  of  subsection  7.3,  was  used  in 
[10,  26]  and  in  chapter  5  of  [6].  Similar  to  subsection  7.3,  only  one  peak  for  each  time-resolved 
curve  was  singled  out  in  [6,  10,  26].  In  that  case  transmitted  (rather  than  backscattering)  time- 
resolved  experimental  data  were  measured  on  many  detectors  that  were  located  on  a  piece  of 
a  plane.  Because  those  data  were  collected  on  a  piece  of  a  plane,  it  was  possible  to  compute 
3D  images  by  a  3D  analogue  of  the  algorithm  of  this  paper.  The  data  of  [6,  10,  26]  were 
collected  in  a  controlled  laboratory  environment,  unlike  the  data  of  this  paper,  which  were 
collected  in  the  uncontrolled  environment  in  the  field.  Another  important  element,  which  led 
to  some  differences  with  the  method  of  subsection  7.3,  was  that  the  reference  signal  (i.e.  the 
one  for  the  free  space)  was  measured  in  [6,  10,  26].  The  latter  is  not  our  case.  Also,  unlike 
this  paper,  refractive  indices  n  =  of  targets  were  directly  measured  a  posteriori,  i.e.  after 
computational  results  were  obtained.  In  other  words,  refractive  indices  were  computed  for 
the  most  difficult  blind  data  case.  Very  accurate  reconstructions  of  both  refractive  indices  and 
locations  of  targets  were  obtained  in  [6,  26] .  Tables  5  and  6  in  [26]  as  well  as  tables  5.4  and  5.5 
in  [6]  show  only  a  few  per  cent  difference  between  computed  and  directly  measured  refractive 
indices  in  six  out  of  six  cases.  Next,  the  so-called  adaptivity  technique  [3,  4]  has  refined 
images.  The  adaptivity  has  provided  quite  accurate  reconstructions  of  all  three  components 
of  targets:  refractive  indices,  locations  and  shapes;  see  [10]  and  figures  5.13-5.16  in  [6],  We 
point  out  that  the  adaptivity  has  used  the  solution  of  the  approximately  globally  convergent 
algorithm  as  the  starting  point.  There  is  no  reason  to  apply  the  adaptivity  to  the  data  of  this 
paper,  because  of  the  3D  reality  versus  the  ID  mathematical  model  (see  item  3  above). 

Since  the  data  pre-processing  procedure  of  subsection  7.3  results  only  in  a  single  peak  to 
work  with,  it  might  miss  some  important  yet  unknown  pieces  of  information.  Still,  because 
of  the  analogy  of  this  procedure  with  that  of  [6,  10,  26],  we  believe  that  the  resulting  pre- 
processed  data  of  figure  5(b)  are  sufficient  for  the  very  limited  goal  given  above.  We  believe 
that  table  1  indicates  that  we  have  achieved  that  goal.  It  remains  to  be  seen  whether  the  number 
R  can  be  satisfactory  computed  by  some  other  techniques  if  using  our  data. 
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The  scaling  number  SN  =  1 0  7 ,  which  we  chose  on  the  basis  of  considerations  of 
subsection  7.3,  might  not  be  an  optimal  one.  Hence,  we  conducted  a  limited  sensitivity  study 
to  see  how  sensitive  R  is  to  the  choice  of  scaling  number;  see  test  2  in  section  8.  We  found 
that  variations  of  SN  by  20%,  SN  =  0.8  x  10~7  and  SN  —  1.2  x  10-7  result  in  values  of  R, 
which  are  still  within  tabulated  limits  for  all  five  targets  we  worked  with.  An  optimal  value 
of  SN  might  likely  be  computed  via  comparison  of  values  of  R  R  (SN)  with  the  directly 
measured  value  of  R  for  a  few  targets.  Next,  the  so-chosen  SN  should  be  used  for  all  other 
targets.  At  this  point,  however,  we  are  not  in  a  position  to  do  this  because  of  the  absence  of 
direct  measurements  of  dielectric  constants  of  the  above  targets. 

A  more  complete  understanding  of  the  capability  of  our  algorithm  to  work  with  more 
complicated  targets  as  well  as  a  better  understanding  of  its  accuracy  limits  for  this  kind  of 
experimental  data  will  be  one  of  the  topics  of  our  future  research  effort.  In  particular,  that 
effort  will  likely  include  a  more  sophisticated  data  pre-processing  procedure,  which  would 
result  in  a  bigger  information  content  being  extracted  from  each  experimental  curve. 

The  recovered  dielectric  constant  by  itself  is  not  sufficient  information  to  distinguish  one 
target  from  another.  The  purpose  of  estimating  the  dielectric  constant  is  to  provide  one  extra 
piece  of  information  about  the  target.  Indeed,  up  to  this  point,  most  of  the  radar  community 
relies  solely  on  the  intensity  of  the  radar  image  for  the  detection  and  discrimination  of  targets. 
It  is  hoped  therefore  that  when  the  intensity  information  is  coupled  with  the  new  dielectric 
information,  algorithms  can  then  be  designed  that  will  ultimately  provide  better  performance 
in  terms  of  probability  of  detection  and  false  alarm  rate.  As  is  clear  from  table  1,  some  targets 
will  have  dielectric  values  that  tend  to  group  together,  but  even  that  is  useful  information. 
For  example,  if  the  estimated  dielectric  value  is  consistent  with  a  plastic  land  mine,  then  this 
would  be  another  clue  to  uncovering  the  target. 

In  summary,  we  believe  that  the  results  of  this  paper  indicate  that  our  algorithm  is 
capable  of  stably  and  reliably  calculating  both  the  function  R(x)  and  the  number  R  from  the 
experimental  data  above.  On  the  next  step  of  this  research  one  should  directly  measure  values 
of  dielectric  constants  of  targets  and  then  compare  measured  and  computational  results.  Such 
a  study  might  clarify  at  least  some  questions  discussed  above  in  this  section.  However,  the 
process  of  collection  of  the  above  experimental  data  is  both  expensive  and  time  consuming. 
This  is  the  reason  why  we  do  not  have  any  other  experimental  data  at  the  moment.  Still,  we 
might  well  get  them  in  the  future. 
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